Fraxinus  16.5.0-fx-rc1
An IGT application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
cxDicomConverter.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 "cxDicomConverter.h"
34 #include "ctkDICOMDatabase.h"
35 #include "cxTypeConversions.h"
36 #include "cxTime.h"
37 #include "vtkImageData.h"
39 #include <vtkImageAppend.h>
40 #include <vtkImageCast.h>
41 #include "cxReporter.h"
42 
43 #include "cxLogger.h"
44 #include "ctkDICOMItem.h"
45 #include "dcfilefo.h" // DcmFileFormat
46 #include "dcdeftag.h" // defines all dcm tags
47 #include "dcmimage.h"
48 #include <string.h>
49 
50 #include "cxDicomImageReader.h"
51 
52 typedef vtkSmartPointer<vtkImageAppend> vtkImageAppendPtr;
53 
54 namespace cx
55 {
56 
58 {
59 }
60 
62 {
63 }
64 
65 void DicomConverter::setDicomDatabase(ctkDICOMDatabase* database)
66 {
67  mDatabase = database;
68 }
69 
70 QString DicomConverter::generateUid(DicomImageReaderPtr reader)
71 {
72  QString seriesDescription = reader->item()->GetElementAsString(DCM_SeriesDescription);
73  QString seriesNumber = reader->item()->GetElementAsString(DCM_SeriesNumber);
74 
75  // uid: uid _ <timestamp>
76  // name: find something from series
77  QString currentTimestamp = QDateTime::currentDateTime().toString(timestampSecondsFormat());
78  QString uid = QString("%1_%2_%3").arg(seriesDescription).arg(seriesNumber).arg(currentTimestamp);
79  uid = this->convertToValidFilename(uid);
80  return uid;
81 }
82 
83 QString DicomConverter::convertToValidFilename(QString text) const
84 {
85  QStringList illegal;
86  illegal << "\\s" << "\\." << ":" << ";" << "\\<" << "\\>" << "\\*" << "\\^" << ",";
87  QRegExp regexp(QString("(%1)").arg(illegal.join("|")));
88  text = text.replace(regexp, "_");
89  return text ;
90 }
91 
92 QString DicomConverter::generateName(DicomImageReaderPtr reader)
93 {
94  QString seriesDescription = reader->item()->GetElementAsString(DCM_SeriesDescription);
95  QString name = QString("%1").arg(seriesDescription);
96  return name;
97 }
98 
99 ImagePtr DicomConverter::createCxImageFromDicomFile(QString filename, bool ignoreLocalizerImages)
100 {
102  if (!reader)
103  {
104  reportWarning(QString("File not found: %1").arg(filename));
105  return ImagePtr();
106  }
107 
108  if(ignoreLocalizerImages && reader->isLocalizerImage())
109  {
110  reportWarning(QString("Localizer image removed from series: %1").arg(filename));
111  return ImagePtr();
112  }
113 
114  if (reader->getNumberOfFrames()==0)
115  {
116  reportWarning(QString("Found no images in %1, skipping.").arg(filename));
117  return ImagePtr();
118  }
119 
120  QString uid = this->generateUid(reader);
121  QString name = this->generateName(reader);
122  cx::ImagePtr image = cx::Image::create(uid, name);
123 
124  vtkImageDataPtr imageData = reader->createVtkImageData();
125  if (!imageData)
126  {
127  reportWarning(QString("Failed to create image for %1.").arg(filename));
128  return ImagePtr();
129  }
130  image->setVtkImageData(imageData);
131 
132  QString modality = reader->item()->GetElementAsString(DCM_Modality);
133  image->setModality(modality);
134 
135  DicomImageReader::WindowLevel windowLevel = reader->getWindowLevel();
136  image->setInitialWindowLevel(windowLevel.width, windowLevel.center);
137 
138  Transform3D M = reader->getImageTransformPatient();
139  image->get_rMd_History()->setRegistration(M);
140 
141  reportDebug(QString("Image created from %1").arg(filename));
142  return image;
143 }
144 
145 std::vector<ImagePtr> DicomConverter::createImages(QStringList files)
146 {
147  std::vector<ImagePtr> retval;
148  for (int i=0; i<files.size(); ++i)
149  {
150  bool ignoreSpesialImages = true;
151  ImagePtr image = this->createCxImageFromDicomFile(files[i], ignoreSpesialImages);
152  if (image)
153  retval.push_back(image);
154  }
155  return retval;
156 }
157 
158 std::map<double, ImagePtr> DicomConverter::sortImagesAlongDirection(std::vector<ImagePtr> images, Vector3D e_sort)
159 {
160  std::map<double, ImagePtr> sorted;
161  for (int i=0; i<images.size(); ++i)
162  {
163  Vector3D pos = images[i]->get_rMd().coord(Vector3D(0,0,0));
164  double dist = dot(pos, e_sort);
165 
166  sorted[dist] = images[i];
167  }
168  return sorted;
169 }
170 
171 bool DicomConverter::slicesFormRegularGrid(std::map<double, ImagePtr> sorted, Vector3D e_sort) const
172 {
173  std::vector<Vector3D> positions;
174  std::vector<double> distances;
175  for (std::map<double, ImagePtr>::iterator iter=sorted.begin(); iter!=sorted.end(); ++iter)
176  {
177  ImagePtr current = iter->second;
178 
179  Vector3D pos = current->get_rMd().coord(Vector3D(0,0,0));
180  positions.push_back(pos);
181 
182  if (positions.size()>=2)
183  {
184  Vector3D p0 = positions[positions.size()-2];
185  Vector3D p1 = positions[positions.size()-1];
186  double dist = dot(p1-p0, e_sort);
187  distances.push_back(dist);
188 
189  Vector3D tilt = cross(p1-p0, e_sort);
190  double sliceGantryTiltTolerance = 0.001;
191  if (!similar(tilt.length(), 0.0, sliceGantryTiltTolerance))
192  {
193  reportError(QString("Dicom convert: found gantry tilt: %1, cannot create image.").arg(tilt.length()));
194  return false;
195  }
196  }
197 
198  if (distances.size()>=2)
199  {
200  double d0 = distances[distances.size()-2];
201  double d1 = distances[distances.size()-1];
202  double sliceSpacingTolerance = 0.01;
203  if (!similar(d0, d1, sliceSpacingTolerance))
204  {
205  reportError(QString("Dicom convert: found uneven slice spacing: %1 != %2, cannot create image.").arg(d0).arg(d1));
206  return false;
207  }
208  }
209  }
210 
211  return true;
212 }
213 
214 double DicomConverter::getMeanSliceDistance(std::map<double, ImagePtr> sorted) const
215 {
216  if (sorted.size()==0)
217  return 0;
218 
219  // check for multislice image
220  vtkImageDataPtr first = sorted.begin()->second->getBaseVtkImageData();
221  if (first->GetDimensions()[2]>1)
222  return first->GetSpacing()[2];
223 
224  if (sorted.size()<2)
225  return 0;
226 
227  // use average of all slices
228  double p1 = sorted.rbegin()->first;
229  double p0 = sorted.begin()->first;
230  return (p1-p0)/sorted.size();
231 }
232 
233 ImagePtr DicomConverter::mergeSlices(std::map<double, ImagePtr> sorted) const
234 {
235  vtkImageAppendPtr appender = vtkImageAppendPtr::New();
236  appender->SetAppendAxis(2);
237 
238  ImagePtr retval = sorted.begin()->second;
239 
240  int i = 0;
241 
242  for (std::map<double, ImagePtr>::iterator iter=sorted.begin(); iter!=sorted.end(); ++iter)
243  {
244  ImagePtr current = iter->second;
245 
246  // Set window width and level to the values of the middle frame
247  if (i == sorted.size() / 2)
248  retval->setInitialWindowLevel(current->getInitialWindowWidth(), current->getInitialWindowLevel());
249  ++i;
250 
251  //Convert all slices to same format
252  vtkImageCastPtr imageCast = vtkImageCastPtr::New();
253  imageCast->SetInputData(current->getBaseVtkImageData());
254  imageCast->SetOutputScalarTypeToShort();
255  imageCast->Update();
256 
257  appender->AddInputData(imageCast->GetOutput());
258  }
259  appender->Update();
260 
261  vtkImageDataPtr wholeImage = appender->GetOutput();
262  Eigen::Array3d spacing(wholeImage->GetSpacing());
263  spacing[2] = this->getMeanSliceDistance(sorted);
264  wholeImage->SetSpacing(spacing.data());
265 
266  retval->setVtkImageData(wholeImage);
267 
268  return retval;
269 }
270 
272 {
273  QStringList files = mDatabase->filesForSeries(series);
274 
275  std::vector<ImagePtr> images = this->createImages(files);
276 
277  if (images.empty())
278  return ImagePtr();
279 
280  if (images.size()==1)
281  {
282  return images.front();
283  }
284 
285  Vector3D e_sort = images.front()->get_rMd().vector(Vector3D(0,0,1));
286 
287  std::map<double, ImagePtr> sorted = this->sortImagesAlongDirection(images, e_sort);
288 
289  if (!this->slicesFormRegularGrid(sorted, e_sort))
290  return ImagePtr();
291 
292  ImagePtr retval = this->mergeSlices(sorted);
293  return retval;
294 }
295 
296 } /* namespace cx */
static DicomImageReaderPtr createFromFile(QString filename)
void reportError(QString msg)
Definition: cxLogger.cpp:92
static ImagePtr create(const QString &uid, const QString &name)
Definition: cxImage.cpp:118
Transform3D Transform3D
Transform3D is a representation of an affine 3D transform.
boost::shared_ptr< class Image > ImagePtr
Definition: cxDicomWidget.h:48
void setDicomDatabase(ctkDICOMDatabase *database)
QString timestampSecondsFormat()
Definition: cxTime.cpp:39
bool similar(const DoubleBoundingBox3D &a, const DoubleBoundingBox3D &b, double tol)
Vector3D cross(const Vector3D &a, const Vector3D &b)
compute cross product of a and b.
Definition: cxVector3D.cpp:62
ImagePtr convertToImage(QString seriesUid)
void reportWarning(QString msg)
Definition: cxLogger.cpp:91
double dot(const Vector3D &a, const Vector3D &b)
compute inner product (or dot product) of a and b.
Definition: cxVector3D.cpp:67
vtkSmartPointer< class vtkImageCast > vtkImageCastPtr
boost::shared_ptr< class DicomImageReader > DicomImageReaderPtr
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
Definition: cxVector3D.h:63
vtkSmartPointer< class vtkImageData > vtkImageDataPtr
void reportDebug(QString msg)
Definition: cxLogger.cpp:89
vtkSmartPointer< vtkImageAppend > vtkImageAppendPtr