34 #include "ctkDICOMDatabase.h"
37 #include "vtkImageData.h"
39 #include <vtkImageAppend.h>
40 #include <vtkImageCast.h>
44 #include "ctkDICOMItem.h"
72 QString seriesDescription = reader->item()->GetElementAsString(DCM_SeriesDescription);
73 QString seriesNumber = reader->item()->GetElementAsString(DCM_SeriesNumber);
78 QString uid = QString(
"%1_%2_%3").arg(seriesDescription).arg(seriesNumber).arg(currentTimestamp);
79 uid = this->convertToValidFilename(uid);
83 QString DicomConverter::convertToValidFilename(QString text)
const
86 illegal <<
"\\s" <<
"\\." <<
":" <<
";" <<
"\\<" <<
"\\>" <<
"\\*" <<
"\\^" <<
",";
87 QRegExp regexp(QString(
"(%1)").arg(illegal.join(
"|")));
88 text = text.replace(regexp,
"_");
94 QString seriesDescription = reader->item()->GetElementAsString(DCM_SeriesDescription);
95 QString name = QString(
"%1").arg(seriesDescription);
99 ImagePtr DicomConverter::createCxImageFromDicomFile(QString filename,
bool ignoreLocalizerImages)
108 if(ignoreLocalizerImages && reader->isLocalizerImage())
110 reportWarning(QString(
"Localizer image removed from series: %1").arg(filename));
114 if (reader->getNumberOfFrames()==0)
116 reportWarning(QString(
"Found no images in %1, skipping.").arg(filename));
120 QString uid = this->generateUid(reader);
121 QString name = this->generateName(reader);
127 reportWarning(QString(
"Failed to create image for %1.").arg(filename));
130 image->setVtkImageData(imageData);
132 QString modality = reader->item()->GetElementAsString(DCM_Modality);
133 image->setModality(modality);
135 DicomImageReader::WindowLevel windowLevel = reader->getWindowLevel();
136 image->setInitialWindowLevel(windowLevel.width, windowLevel.center);
138 Transform3D M = reader->getImageTransformPatient();
139 image->get_rMd_History()->setRegistration(M);
141 reportDebug(QString(
"Image created from %1").arg(filename));
145 std::vector<ImagePtr> DicomConverter::createImages(QStringList files)
147 std::vector<ImagePtr> retval;
148 for (
int i=0; i<files.size(); ++i)
150 bool ignoreSpesialImages =
true;
151 ImagePtr image = this->createCxImageFromDicomFile(files[i], ignoreSpesialImages);
153 retval.push_back(image);
158 std::map<double, ImagePtr> DicomConverter::sortImagesAlongDirection(std::vector<ImagePtr> images,
Vector3D e_sort)
160 std::map<double, ImagePtr> sorted;
161 for (
int i=0; i<images.size(); ++i)
164 double dist =
dot(pos, e_sort);
166 sorted[dist] = images[i];
171 bool DicomConverter::slicesFormRegularGrid(std::map<double, ImagePtr> sorted,
Vector3D e_sort)
const
173 std::vector<Vector3D> positions;
174 std::vector<double> distances;
175 for (std::map<double, ImagePtr>::iterator iter=sorted.begin(); iter!=sorted.end(); ++iter)
180 positions.push_back(pos);
182 if (positions.size()>=2)
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);
190 double sliceGantryTiltTolerance = 0.001;
191 if (!
similar(tilt.length(), 0.0, sliceGantryTiltTolerance))
193 reportError(QString(
"Dicom convert: found gantry tilt: %1, cannot create image.").arg(tilt.length()));
198 if (distances.size()>=2)
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))
205 reportError(QString(
"Dicom convert: found uneven slice spacing: %1 != %2, cannot create image.").arg(d0).arg(d1));
214 double DicomConverter::getMeanSliceDistance(std::map<double, ImagePtr> sorted)
const
216 if (sorted.size()==0)
221 if (first->GetDimensions()[2]>1)
222 return first->GetSpacing()[2];
228 double p1 = sorted.rbegin()->first;
229 double p0 = sorted.begin()->first;
230 return (p1-p0)/sorted.size();
233 ImagePtr DicomConverter::mergeSlices(std::map<double, ImagePtr> sorted)
const
236 appender->SetAppendAxis(2);
238 ImagePtr retval = sorted.begin()->second;
242 for (std::map<double, ImagePtr>::iterator iter=sorted.begin(); iter!=sorted.end(); ++iter)
247 if (i == sorted.size() / 2)
248 retval->setInitialWindowLevel(current->getInitialWindowWidth(), current->getInitialWindowLevel());
253 imageCast->SetInputData(current->getBaseVtkImageData());
254 imageCast->SetOutputScalarTypeToShort();
257 appender->AddInputData(imageCast->GetOutput());
262 Eigen::Array3d spacing(wholeImage->GetSpacing());
263 spacing[2] = this->getMeanSliceDistance(sorted);
264 wholeImage->SetSpacing(spacing.data());
266 retval->setVtkImageData(wholeImage);
273 QStringList files = mDatabase->filesForSeries(series);
275 std::vector<ImagePtr> images = this->createImages(files);
280 if (images.size()==1)
282 return images.front();
287 std::map<double, ImagePtr> sorted = this->sortImagesAlongDirection(images, e_sort);
289 if (!this->slicesFormRegularGrid(sorted, e_sort))
292 ImagePtr retval = this->mergeSlices(sorted);
static DicomImageReaderPtr createFromFile(QString filename)
void reportError(QString msg)
static ImagePtr create(const QString &uid, const QString &name)
Transform3D Transform3D
Transform3D is a representation of an affine 3D transform.
boost::shared_ptr< class Image > ImagePtr
void setDicomDatabase(ctkDICOMDatabase *database)
QString timestampSecondsFormat()
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.
ImagePtr convertToImage(QString seriesUid)
void reportWarning(QString msg)
virtual ~DicomConverter()
double dot(const Vector3D &a, const Vector3D &b)
compute inner product (or dot product) of a and b.
vtkSmartPointer< class vtkImageCast > vtkImageCastPtr
boost::shared_ptr< class DicomImageReader > DicomImageReaderPtr
Eigen::Vector3d Vector3D
Vector3D is a representation of a point or vector in 3D.
void reportDebug(QString msg)
vtkSmartPointer< vtkImageAppend > vtkImageAppendPtr