diff --git a/Core/Code/IO/mitkDicomSR_GantryTiltInformation.cpp b/Core/Code/IO/mitkDicomSR_GantryTiltInformation.cpp new file mode 100644 index 0000000000..6a324d2b64 --- /dev/null +++ b/Core/Code/IO/mitkDicomSR_GantryTiltInformation.cpp @@ -0,0 +1,201 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include + +namespace mitk +{ + +DicomSeriesReader::GantryTiltInformation::GantryTiltInformation() +: m_ShiftUp(0.0) +, m_ShiftRight(0.0) +, m_ShiftNormal(0.0) +, m_ITKAssumedSliceSpacing(0.0) +, m_NumberOfSlicesApart(1) +{ +} + + +#define doublepoint(x) \ + Point3Dd x; \ + x[0] = x ## f[0]; \ + x[1] = x ## f[1]; \ + x[2] = x ## f[2]; + + +#define doublevector(x) \ + Vector3Dd x; \ + x[0] = x ## f[0]; \ + x[1] = x ## f[1]; \ + x[2] = x ## f[2]; + +DicomSeriesReader::GantryTiltInformation::GantryTiltInformation( + const Point3D& origin1f, const Point3D& origin2f, + const Vector3D& rightf, const Vector3D& upf, + unsigned int numberOfSlicesApart) +: m_ShiftUp(0.0) +, m_ShiftRight(0.0) +, m_ShiftNormal(0.0) +, m_NumberOfSlicesApart(numberOfSlicesApart) +{ + assert(numberOfSlicesApart); + + doublepoint(origin1); + doublepoint(origin2); + doublevector(right); + doublevector(up); + + // determine if slice 1 (imagePosition1 and imageOrientation1) and slice 2 can be in one orthogonal slice stack: + // calculate a line from origin 1, directed along the normal of slice (calculated as the cross product of orientation 1) + // check if this line passes through origin 2 + + /* + Determine if line (imagePosition2 + l * normal) contains imagePosition1. + Done by calculating the distance of imagePosition1 from line (imagePosition2 + l *normal) + + E.g. http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html + + squared distance = | (pointAlongNormal - origin2) x (origin2 - origin1) | ^ 2 + / + |pointAlongNormal - origin2| ^ 2 + + ( x meaning the cross product ) + */ + + Vector3Dd normal = itk::CrossProduct(right, up); + Point3Dd pointAlongNormal = origin2 + normal; + + double numerator = itk::CrossProduct( pointAlongNormal - origin2 , origin2 - origin1 ).GetSquaredNorm(); + double denominator = (pointAlongNormal - origin2).GetSquaredNorm(); + + double distance = sqrt(numerator / denominator); + + if ( distance > 0.001 ) // mitk::eps is too small; 1/1000 of a mm should be enough to detect tilt + { + MITK_DEBUG << " Series seems to contain a tilted (or sheared) geometry"; + MITK_DEBUG << " Distance of expected slice origin from actual slice origin: " << distance; + MITK_DEBUG << " ==> storing this shift for later analysis:"; + MITK_DEBUG << " v right: " << right; + MITK_DEBUG << " v up: " << up; + MITK_DEBUG << " v normal: " << normal; + + Point3Dd projectionRight = projectPointOnLine( origin1, origin2, right ); + Point3Dd projectionNormal = projectPointOnLine( origin1, origin2, normal ); + + m_ShiftRight = (projectionRight - origin2).GetNorm(); + m_ShiftNormal = (projectionNormal - origin2).GetNorm(); + + /* + now also check to which side the image is shifted. + + Calculation e.g. from + http://mathworld.wolfram.com/Point-PlaneDistance.html + */ + + Point3Dd testPoint = origin1; + Vector3Dd planeNormal = up; + + double signedDistance = ( + planeNormal[0] * testPoint[0] + + planeNormal[1] * testPoint[1] + + planeNormal[2] * testPoint[2] + - ( + planeNormal[0] * origin2[0] + + planeNormal[1] * origin2[1] + + planeNormal[2] * origin2[2] + ) + ) + / + sqrt( planeNormal[0] * planeNormal[0] + + planeNormal[1] * planeNormal[1] + + planeNormal[2] * planeNormal[2] + ); + + m_ShiftUp = signedDistance; + + m_ITKAssumedSliceSpacing = (origin2 - origin1).GetNorm(); + // How do we now this is assumed? See header documentation for ITK code references + //double itkAssumedSliceSpacing = sqrt( m_ShiftUp * m_ShiftUp + m_ShiftNormal * m_ShiftNormal ); + + MITK_DEBUG << " shift normal: " << m_ShiftNormal; + MITK_DEBUG << " shift normal assumed by ITK: " << m_ITKAssumedSliceSpacing; + MITK_DEBUG << " shift up: " << m_ShiftUp; + MITK_DEBUG << " shift right: " << m_ShiftRight; + + MITK_DEBUG << " tilt angle (deg): " << atan( m_ShiftUp / m_ShiftNormal ) * 180.0 / 3.1415926535; + } +} + +Point3D +DicomSeriesReader::GantryTiltInformation::projectPointOnLine( Point3Dd p, Point3Dd lineOrigin, Vector3Dd lineDirection ) +{ + /** + See illustration at http://mo.mathematik.uni-stuttgart.de/inhalt/aussage/aussage472/ + + vector(lineOrigin,p) = normal * ( innerproduct((p - lineOrigin),normal) / squared-length(normal) ) + */ + + Vector3Dd lineOriginToP = p - lineOrigin; + double innerProduct = lineOriginToP * lineDirection; + + double factor = innerProduct / lineDirection.GetSquaredNorm(); + Point3Dd projection = lineOrigin + factor * lineDirection; + + return projection; +} + +double +DicomSeriesReader::GantryTiltInformation::GetTiltCorrectedAdditionalSize() const +{ + return fabs(m_ShiftUp); +} + +double +DicomSeriesReader::GantryTiltInformation::GetTiltAngleInDegrees() const +{ + return atan( fabs(m_ShiftUp) / m_ShiftNormal ) * 180.0 / 3.1415926535; +} + +double +DicomSeriesReader::GantryTiltInformation::GetMatrixCoefficientForCorrectionInWorldCoordinates() const +{ + // so many mm need to be shifted per slice! + return m_ShiftUp / static_cast(m_NumberOfSlicesApart); +} + +double +DicomSeriesReader::GantryTiltInformation::GetRealZSpacing() const +{ + return m_ShiftNormal / static_cast(m_NumberOfSlicesApart); +} + + +bool +DicomSeriesReader::GantryTiltInformation::IsSheared() const +{ + return ( fabs(m_ShiftRight) > 0.001 + || fabs(m_ShiftUp) > 0.001); +} + + +bool +DicomSeriesReader::GantryTiltInformation::IsRegularGantryTilt() const +{ + return ( fabs(m_ShiftRight) < 0.001 + && fabs(m_ShiftUp) > 0.001); +} + +} // end namespace mitk diff --git a/Core/Code/IO/mitkDicomSR_ImageBlockDescriptor.cpp b/Core/Code/IO/mitkDicomSR_ImageBlockDescriptor.cpp new file mode 100644 index 0000000000..ce637f8192 --- /dev/null +++ b/Core/Code/IO/mitkDicomSR_ImageBlockDescriptor.cpp @@ -0,0 +1,242 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include + +#include + +namespace mitk +{ + +DicomSeriesReader::ImageBlockDescriptor::ImageBlockDescriptor() +:m_HasGantryTiltCorrected(false) +,m_HasMultipleTimePoints(false) +,m_IsMultiFrameImage(false) +{ +} + +DicomSeriesReader::ImageBlockDescriptor::~ImageBlockDescriptor() +{ + // nothing +} + +DicomSeriesReader::ImageBlockDescriptor::ImageBlockDescriptor(const StringContainer& files) +:m_HasGantryTiltCorrected(false) +,m_HasMultipleTimePoints(false) +,m_IsMultiFrameImage(false) +{ + m_Filenames = files; +} + + +void DicomSeriesReader::ImageBlockDescriptor::AddFile(const std::string& filename) +{ + m_Filenames.push_back( filename ); +} + +void DicomSeriesReader::ImageBlockDescriptor::AddFiles(const StringContainer& files) +{ + m_Filenames.insert( m_Filenames.end(), files.begin(), files.end() ); +} + +DicomSeriesReader::StringContainer DicomSeriesReader::ImageBlockDescriptor::GetFilenames() const +{ + return m_Filenames; +} + +std::string DicomSeriesReader::ImageBlockDescriptor::GetImageBlockUID() const +{ + return m_ImageBlockUID; +} + +std::string DicomSeriesReader::ImageBlockDescriptor::GetSeriesInstanceUID() const +{ + return m_SeriesInstanceUID; +} + +std::string DicomSeriesReader::ImageBlockDescriptor::GetModality() const +{ + return m_Modality; +} + +std::string DicomSeriesReader::ImageBlockDescriptor::GetSOPClassUIDAsString() const +{ + gdcm::UIDs uidKnowledge; + uidKnowledge.SetFromUID( m_SOPClassUID.c_str() ); + return uidKnowledge.GetName(); +} + +std::string DicomSeriesReader::ImageBlockDescriptor::GetSOPClassUID() const +{ + return m_SOPClassUID; +} + +bool DicomSeriesReader::ImageBlockDescriptor::IsMultiFrameImage() const +{ + return m_IsMultiFrameImage; +} + +DicomSeriesReader::ReaderImplementationLevel DicomSeriesReader::ImageBlockDescriptor::GetReaderImplementationLevel() const +{ + if ( this->IsMultiFrameImage() ) + return ReaderImplementationLevel_Unsupported; + + gdcm::UIDs uidKnowledge; + uidKnowledge.SetFromUID( m_SOPClassUID.c_str() ); + + gdcm::UIDs::TSType uid = uidKnowledge; + + switch (uid) + { + case gdcm::UIDs::CTImageStorage: + case gdcm::UIDs::MRImageStorage: + case gdcm::UIDs::PositronEmissionTomographyImageStorage: + case gdcm::UIDs::ComputedRadiographyImageStorage: + case gdcm::UIDs::DigitalXRayImageStorageForPresentation: + case gdcm::UIDs::DigitalXRayImageStorageForProcessing: + return ReaderImplementationLevel_Supported; + + case gdcm::UIDs::NuclearMedicineImageStorage: + return ReaderImplementationLevel_PartlySupported; + + case gdcm::UIDs::SecondaryCaptureImageStorage: + return ReaderImplementationLevel_Implemented; + + default: + return ReaderImplementationLevel_Unsupported; + } +} + +bool DicomSeriesReader::ImageBlockDescriptor::HasGantryTiltCorrected() const +{ + return m_HasGantryTiltCorrected; +} + +/* + PS defined IPS defined PS==IPS + 0 0 --> UNKNOWN spacing, loader will invent + 0 1 --> spacing as at detector surface + 1 0 --> spacing as in patient + 1 1 0 --> detector surface spacing CORRECTED for geometrical magnifications: spacing as in patient + 1 1 1 --> detector surface spacing NOT corrected for geometrical magnifications: spacing as at detector +*/ +DicomSeriesReader::PixelSpacingInterpretation DicomSeriesReader::ImageBlockDescriptor::GetPixelSpacingType() const +{ + if (m_PixelSpacing.empty()) + { + if (m_ImagerPixelSpacing.empty()) + { + return PixelSpacingInterpretation_SpacingUnknown; + } + else + { + return PixelSpacingInterpretation_SpacingAtDetector; + } + } + else // Pixel Spacing defined + { + if (m_ImagerPixelSpacing.empty()) + { + return PixelSpacingInterpretation_SpacingInPatient; + } + else if (m_PixelSpacing != m_ImagerPixelSpacing) + { + return PixelSpacingInterpretation_SpacingInPatient; + } + else + { + return PixelSpacingInterpretation_SpacingAtDetector; + } + } +} + +bool DicomSeriesReader::ImageBlockDescriptor::PixelSpacingRelatesToPatient() const +{ + return GetPixelSpacingType() == PixelSpacingInterpretation_SpacingInPatient; +} + +bool DicomSeriesReader::ImageBlockDescriptor::PixelSpacingRelatesToDetector() const +{ + return GetPixelSpacingType() == PixelSpacingInterpretation_SpacingAtDetector; +} + +bool DicomSeriesReader::ImageBlockDescriptor::PixelSpacingIsUnknown() const +{ + return GetPixelSpacingType() == PixelSpacingInterpretation_SpacingUnknown; +} + +void DicomSeriesReader::ImageBlockDescriptor::SetPixelSpacingInformation(const std::string& pixelSpacing, const std::string& imagerPixelSpacing) +{ + m_PixelSpacing = pixelSpacing; + m_ImagerPixelSpacing = imagerPixelSpacing; +} + +void DicomSeriesReader::ImageBlockDescriptor::GetDesiredMITKImagePixelSpacing( float& spacingX, float& spacingY) const +{ + // preference for "in patient" pixel spacing + if ( !DICOMStringToSpacing( m_PixelSpacing, spacingX, spacingY ) ) + { + // fallback to "on detector" spacing + if ( !DICOMStringToSpacing( m_ImagerPixelSpacing, spacingX, spacingY ) ) + { + // last resort: invent something + spacingX = spacingY = 1.0; + } + } +} + +bool DicomSeriesReader::ImageBlockDescriptor::HasMultipleTimePoints() const +{ + return m_HasMultipleTimePoints; +} + +void DicomSeriesReader::ImageBlockDescriptor::SetImageBlockUID(const std::string& uid) +{ + m_ImageBlockUID = uid; +} + +void DicomSeriesReader::ImageBlockDescriptor::SetSeriesInstanceUID(const std::string& uid) +{ + m_SeriesInstanceUID = uid; +} + +void DicomSeriesReader::ImageBlockDescriptor::SetModality(const std::string& modality) +{ + m_Modality = modality; +} + +void DicomSeriesReader::ImageBlockDescriptor::SetNumberOfFrames(const std::string& numberOfFrames) +{ + m_IsMultiFrameImage = !numberOfFrames.empty(); +} + +void DicomSeriesReader::ImageBlockDescriptor::SetSOPClassUID(const std::string& sopClassUID) +{ + m_SOPClassUID = sopClassUID; +} + + +void DicomSeriesReader::ImageBlockDescriptor::SetHasGantryTiltCorrected(bool on) +{ + m_HasGantryTiltCorrected = on; +} + +void DicomSeriesReader::ImageBlockDescriptor::SetHasMultipleTimePoints(bool on) +{ + m_HasMultipleTimePoints = on; +} + +} // end namespace mitk diff --git a/Core/Code/IO/mitkDicomSR_LoadDICOMRGBPixel.cpp b/Core/Code/IO/mitkDicomSR_LoadDICOMRGBPixel.cpp new file mode 100644 index 0000000000..55e4e38ddd --- /dev/null +++ b/Core/Code/IO/mitkDicomSR_LoadDICOMRGBPixel.cpp @@ -0,0 +1,54 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "mitkDicomSeriesReader.txx" + +namespace mitk +{ + +Image::Pointer +DicomSeriesReader +::MultiplexLoadDICOMByITKRGBPixel(const StringContainer& filenames, bool correctTilt, const GantryTiltInformation& tiltInfo, DcmIoType::Pointer& io, CallbackCommand* command, Image::Pointer preLoadedImageBlock) +{ + switch (io->GetComponentType()) + { + case DcmIoType::UCHAR: + return LoadDICOMByITK< itk::RGBPixel >(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::CHAR: + return LoadDICOMByITK >(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::USHORT: + return LoadDICOMByITK >(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::SHORT: + return LoadDICOMByITK >(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::UINT: + return LoadDICOMByITK >(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::INT: + return LoadDICOMByITK >(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::ULONG: + return LoadDICOMByITK >(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::LONG: + return LoadDICOMByITK >(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::FLOAT: + return LoadDICOMByITK >(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::DOUBLE: + return LoadDICOMByITK >(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + default: + MITK_ERROR << "Found unsupported DICOM scalar pixel type: (enum value) " << io->GetComponentType(); + return NULL; + } +} + +} // end namespace mitk diff --git a/Core/Code/IO/mitkDicomSR_LoadDICOMRGBPixel4D.cpp b/Core/Code/IO/mitkDicomSR_LoadDICOMRGBPixel4D.cpp new file mode 100644 index 0000000000..caad59e3e2 --- /dev/null +++ b/Core/Code/IO/mitkDicomSR_LoadDICOMRGBPixel4D.cpp @@ -0,0 +1,54 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "mitkDicomSeriesReader.txx" + +namespace mitk +{ + +Image::Pointer +DicomSeriesReader +::MultiplexLoadDICOMByITK4DRGBPixel( std::list& imageBlocks, ImageBlockDescriptor imageBlockDescriptor, bool correctTilt, const GantryTiltInformation& tiltInfo, DcmIoType::Pointer& io, CallbackCommand* command, Image::Pointer preLoadedImageBlock) +{ + switch (io->GetComponentType()) + { + case DcmIoType::UCHAR: + return LoadDICOMByITK4D< itk::RGBPixel >(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::CHAR: + return LoadDICOMByITK4D >(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::USHORT: + return LoadDICOMByITK4D >(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::SHORT: + return LoadDICOMByITK4D >(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::UINT: + return LoadDICOMByITK4D >(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::INT: + return LoadDICOMByITK4D >(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::ULONG: + return LoadDICOMByITK4D >(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::LONG: + return LoadDICOMByITK4D >(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::FLOAT: + return LoadDICOMByITK4D >(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::DOUBLE: + return LoadDICOMByITK4D >(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + default: + MITK_ERROR << "Found unsupported DICOM scalar pixel type: (enum value) " << io->GetComponentType(); + return NULL; + } +} + +} // end namespace mitk diff --git a/Core/Code/IO/mitkDicomSR_LoadDICOMScalar.cpp b/Core/Code/IO/mitkDicomSR_LoadDICOMScalar.cpp new file mode 100644 index 0000000000..b469205b30 --- /dev/null +++ b/Core/Code/IO/mitkDicomSR_LoadDICOMScalar.cpp @@ -0,0 +1,57 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "mitkDicomSeriesReader.txx" + +namespace mitk +{ + +Image::Pointer +DicomSeriesReader +::MultiplexLoadDICOMByITKScalar(const StringContainer& filenames, bool correctTilt, const GantryTiltInformation& tiltInfo, DcmIoType::Pointer& io, CallbackCommand* command, Image::Pointer preLoadedImageBlock) +{ + switch (io->GetComponentType()) + { + case DcmIoType::UCHAR: + return LoadDICOMByITK(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::CHAR: + return LoadDICOMByITK(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::USHORT: + return LoadDICOMByITK(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::SHORT: + return LoadDICOMByITK(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::UINT: + return LoadDICOMByITK(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::INT: + return LoadDICOMByITK(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::ULONG: + return LoadDICOMByITK(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::LONG: + return LoadDICOMByITK(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::FLOAT: + return LoadDICOMByITK(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::DOUBLE: + return LoadDICOMByITK(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + default: + MITK_ERROR << "Found unsupported DICOM scalar pixel type: (enum value) " << io->GetComponentType(); + return NULL; + } +} + +} // end namespace mitk + +#include + diff --git a/Core/Code/IO/mitkDicomSR_LoadDICOMScalar4D.cpp b/Core/Code/IO/mitkDicomSR_LoadDICOMScalar4D.cpp new file mode 100644 index 0000000000..953d9a32d0 --- /dev/null +++ b/Core/Code/IO/mitkDicomSR_LoadDICOMScalar4D.cpp @@ -0,0 +1,57 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include "mitkDicomSeriesReader.txx" + +namespace mitk +{ + +Image::Pointer +DicomSeriesReader +::MultiplexLoadDICOMByITK4DScalar( std::list& imageBlocks, ImageBlockDescriptor imageBlockDescriptor, bool correctTilt, const GantryTiltInformation& tiltInfo, DcmIoType::Pointer& io, CallbackCommand* command, Image::Pointer preLoadedImageBlock) +{ + switch (io->GetComponentType()) + { + case DcmIoType::UCHAR: + return LoadDICOMByITK4D(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::CHAR: + return LoadDICOMByITK4D(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::USHORT: + return LoadDICOMByITK4D(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::SHORT: + return LoadDICOMByITK4D(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::UINT: + return LoadDICOMByITK4D(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::INT: + return LoadDICOMByITK4D(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::ULONG: + return LoadDICOMByITK4D(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::LONG: + return LoadDICOMByITK4D(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::FLOAT: + return LoadDICOMByITK4D(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + case DcmIoType::DOUBLE: + return LoadDICOMByITK4D(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + default: + MITK_ERROR << "Found unsupported DICOM scalar pixel type: (enum value) " << io->GetComponentType(); + return NULL; + } +} + +} // end namespace mitk + +#include + diff --git a/Core/Code/IO/mitkDicomSR_SliceGroupingResult.cpp b/Core/Code/IO/mitkDicomSR_SliceGroupingResult.cpp new file mode 100644 index 0000000000..e853bb5604 --- /dev/null +++ b/Core/Code/IO/mitkDicomSR_SliceGroupingResult.cpp @@ -0,0 +1,69 @@ +/*=================================================================== + +The Medical Imaging Interaction Toolkit (MITK) + +Copyright (c) German Cancer Research Center, +Division of Medical and Biological Informatics. +All rights reserved. + +This software is distributed WITHOUT ANY WARRANTY; without +even the implied warranty of MERCHANTABILITY or FITNESS FOR +A PARTICULAR PURPOSE. + +See LICENSE.txt or http://www.mitk.org for details. + +===================================================================*/ + +#include + +namespace mitk +{ + +DicomSeriesReader::SliceGroupingAnalysisResult::SliceGroupingAnalysisResult() +:m_GantryTilt(false) +{ +} + +DicomSeriesReader::StringContainer DicomSeriesReader::SliceGroupingAnalysisResult::GetBlockFilenames() +{ + return m_GroupedFiles; +} + +DicomSeriesReader::StringContainer DicomSeriesReader::SliceGroupingAnalysisResult::GetUnsortedFilenames() +{ + return m_UnsortedFiles; +} + +bool DicomSeriesReader::SliceGroupingAnalysisResult::ContainsGantryTilt() +{ + return m_GantryTilt; +} + +void DicomSeriesReader::SliceGroupingAnalysisResult::AddFileToSortedBlock(const std::string& filename) +{ + m_GroupedFiles.push_back( filename ); +} + +void DicomSeriesReader::SliceGroupingAnalysisResult::AddFileToUnsortedBlock(const std::string& filename) +{ + m_UnsortedFiles.push_back( filename ); +} + +void DicomSeriesReader::SliceGroupingAnalysisResult::AddFilesToUnsortedBlock(const StringContainer& filenames) +{ + m_UnsortedFiles.insert( m_UnsortedFiles.end(), filenames.begin(), filenames.end() ); +} + +void DicomSeriesReader::SliceGroupingAnalysisResult::FlagGantryTilt() +{ + m_GantryTilt = true; +} + +void DicomSeriesReader::SliceGroupingAnalysisResult::UndoPrematureGrouping() +{ + assert( !m_GroupedFiles.empty() ); + m_UnsortedFiles.insert( m_UnsortedFiles.begin(), m_GroupedFiles.back() ); + m_GroupedFiles.pop_back(); +} + +} // end namespace mitk diff --git a/Core/Code/IO/mitkDicomSeriesReader.cpp b/Core/Code/IO/mitkDicomSeriesReader.cpp index e2cab0d2c6..1f97e8f467 100644 --- a/Core/Code/IO/mitkDicomSeriesReader.cpp +++ b/Core/Code/IO/mitkDicomSeriesReader.cpp @@ -1,1981 +1,1779 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ // uncomment for learning more about the internal sorting mechanisms //#define MBILOG_ENABLE_DEBUG #include #include #include #include #include #include #include #include #include #include #include "mitkProperties.h" namespace mitk { -typedef itk::GDCMSeriesFileNames DcmFileNamesGeneratorType; - -DicomSeriesReader::ImageBlockDescriptor::ImageBlockDescriptor() -:m_HasGantryTiltCorrected(false) -,m_HasMultipleTimePoints(false) -,m_IsMultiFrameImage(false) -{ -} - -DicomSeriesReader::ImageBlockDescriptor::~ImageBlockDescriptor() -{ - // nothing -} - -DicomSeriesReader::ImageBlockDescriptor::ImageBlockDescriptor(const StringContainer& files) -:m_HasGantryTiltCorrected(false) -,m_HasMultipleTimePoints(false) -,m_IsMultiFrameImage(false) -{ - m_Filenames = files; -} - - -void DicomSeriesReader::ImageBlockDescriptor::AddFile(const std::string& filename) -{ - m_Filenames.push_back( filename ); -} - -void DicomSeriesReader::ImageBlockDescriptor::AddFiles(const StringContainer& files) -{ - m_Filenames.insert( m_Filenames.end(), files.begin(), files.end() ); -} - -DicomSeriesReader::StringContainer DicomSeriesReader::ImageBlockDescriptor::GetFilenames() const -{ - return m_Filenames; -} - -std::string DicomSeriesReader::ImageBlockDescriptor::GetImageBlockUID() const -{ - return m_ImageBlockUID; -} - -std::string DicomSeriesReader::ImageBlockDescriptor::GetSeriesInstanceUID() const -{ - return m_SeriesInstanceUID; -} - -std::string DicomSeriesReader::ImageBlockDescriptor::GetModality() const -{ - return m_Modality; -} - -std::string DicomSeriesReader::ImageBlockDescriptor::GetSOPClassUIDAsString() const -{ - gdcm::UIDs uidKnowledge; - uidKnowledge.SetFromUID( m_SOPClassUID.c_str() ); - return uidKnowledge.GetName(); -} - -std::string DicomSeriesReader::ImageBlockDescriptor::GetSOPClassUID() const -{ - return m_SOPClassUID; -} - -bool DicomSeriesReader::ImageBlockDescriptor::IsMultiFrameImage() const -{ - return m_IsMultiFrameImage; -} - -DicomSeriesReader::ReaderImplementationLevel DicomSeriesReader::ImageBlockDescriptor::GetReaderImplementationLevel() const -{ - if ( this->IsMultiFrameImage() ) - return ReaderImplementationLevel_Unsupported; - - gdcm::UIDs uidKnowledge; - uidKnowledge.SetFromUID( m_SOPClassUID.c_str() ); - - gdcm::UIDs::TSType uid = uidKnowledge; - - switch (uid) - { - case gdcm::UIDs::CTImageStorage: - case gdcm::UIDs::MRImageStorage: - case gdcm::UIDs::PositronEmissionTomographyImageStorage: - case gdcm::UIDs::ComputedRadiographyImageStorage: - case gdcm::UIDs::DigitalXRayImageStorageForPresentation: - case gdcm::UIDs::DigitalXRayImageStorageForProcessing: - return ReaderImplementationLevel_Supported; - - case gdcm::UIDs::NuclearMedicineImageStorage: - return ReaderImplementationLevel_PartlySupported; - - case gdcm::UIDs::SecondaryCaptureImageStorage: - return ReaderImplementationLevel_Implemented; - - default: - return ReaderImplementationLevel_Unsupported; - } -} - std::string DicomSeriesReader::ReaderImplementationLevelToString( const ReaderImplementationLevel& enumValue ) { switch (enumValue) { case ReaderImplementationLevel_Supported: return "Supported"; case ReaderImplementationLevel_PartlySupported: return "PartlySupported"; case ReaderImplementationLevel_Implemented: return "Implemented"; case ReaderImplementationLevel_Unsupported: return "Unsupported"; default: return ""; }; } std::string DicomSeriesReader::PixelSpacingInterpretationToString( const PixelSpacingInterpretation& enumValue ) { switch (enumValue) { case PixelSpacingInterpretation_SpacingInPatient: return "In Patient"; case PixelSpacingInterpretation_SpacingAtDetector: return "At Detector"; case PixelSpacingInterpretation_SpacingUnknown: return "Unknown spacing"; default: return ""; }; } - -bool DicomSeriesReader::ImageBlockDescriptor::HasGantryTiltCorrected() const -{ - return m_HasGantryTiltCorrected; -} - -/* - PS defined IPS defined PS==IPS - 0 0 --> UNKNOWN spacing, loader will invent - 0 1 --> spacing as at detector surface - 1 0 --> spacing as in patient - 1 1 0 --> detector surface spacing CORRECTED for geometrical magnifications: spacing as in patient - 1 1 1 --> detector surface spacing NOT corrected for geometrical magnifications: spacing as at detector -*/ -DicomSeriesReader::PixelSpacingInterpretation DicomSeriesReader::ImageBlockDescriptor::GetPixelSpacingType() const -{ - if (m_PixelSpacing.empty()) - { - if (m_ImagerPixelSpacing.empty()) - { - return PixelSpacingInterpretation_SpacingUnknown; - } - else - { - return PixelSpacingInterpretation_SpacingAtDetector; - } - } - else // Pixel Spacing defined - { - if (m_ImagerPixelSpacing.empty()) - { - return PixelSpacingInterpretation_SpacingInPatient; - } - else if (m_PixelSpacing != m_ImagerPixelSpacing) - { - return PixelSpacingInterpretation_SpacingInPatient; - } - else - { - return PixelSpacingInterpretation_SpacingAtDetector; - } - } -} - -bool DicomSeriesReader::ImageBlockDescriptor::PixelSpacingRelatesToPatient() const -{ - return GetPixelSpacingType() == PixelSpacingInterpretation_SpacingInPatient; -} - -bool DicomSeriesReader::ImageBlockDescriptor::PixelSpacingRelatesToDetector() const -{ - return GetPixelSpacingType() == PixelSpacingInterpretation_SpacingAtDetector; -} - -bool DicomSeriesReader::ImageBlockDescriptor::PixelSpacingIsUnknown() const -{ - return GetPixelSpacingType() == PixelSpacingInterpretation_SpacingUnknown; -} - -void DicomSeriesReader::ImageBlockDescriptor::SetPixelSpacingInformation(const std::string& pixelSpacing, const std::string& imagerPixelSpacing) -{ - m_PixelSpacing = pixelSpacing; - m_ImagerPixelSpacing = imagerPixelSpacing; -} - -void DicomSeriesReader::ImageBlockDescriptor::GetDesiredMITKImagePixelSpacing( float& spacingX, float& spacingY) const -{ - // preference for "in patient" pixel spacing - if ( !DICOMStringToSpacing( m_PixelSpacing, spacingX, spacingY ) ) - { - // fallback to "on detector" spacing - if ( !DICOMStringToSpacing( m_ImagerPixelSpacing, spacingX, spacingY ) ) - { - // last resort: invent something - spacingX = spacingY = 1.0; - } - } -} - -bool DicomSeriesReader::ImageBlockDescriptor::HasMultipleTimePoints() const -{ - return m_HasMultipleTimePoints; -} - -void DicomSeriesReader::ImageBlockDescriptor::SetImageBlockUID(const std::string& uid) -{ - m_ImageBlockUID = uid; -} - -void DicomSeriesReader::ImageBlockDescriptor::SetSeriesInstanceUID(const std::string& uid) -{ - m_SeriesInstanceUID = uid; -} - -void DicomSeriesReader::ImageBlockDescriptor::SetModality(const std::string& modality) -{ - m_Modality = modality; -} - -void DicomSeriesReader::ImageBlockDescriptor::SetNumberOfFrames(const std::string& numberOfFrames) -{ - m_IsMultiFrameImage = !numberOfFrames.empty(); -} - -void DicomSeriesReader::ImageBlockDescriptor::SetSOPClassUID(const std::string& sopClassUID) -{ - m_SOPClassUID = sopClassUID; -} - - -void DicomSeriesReader::ImageBlockDescriptor::SetHasGantryTiltCorrected(bool on) -{ - m_HasGantryTiltCorrected = on; -} - -void DicomSeriesReader::ImageBlockDescriptor::SetHasMultipleTimePoints(bool on) -{ - m_HasMultipleTimePoints = on; -} - - - -DicomSeriesReader::SliceGroupingAnalysisResult::SliceGroupingAnalysisResult() -:m_GantryTilt(false) -{ -} - -DicomSeriesReader::StringContainer DicomSeriesReader::SliceGroupingAnalysisResult::GetBlockFilenames() -{ - return m_GroupedFiles; -} - -DicomSeriesReader::StringContainer DicomSeriesReader::SliceGroupingAnalysisResult::GetUnsortedFilenames() -{ - return m_UnsortedFiles; -} - -bool DicomSeriesReader::SliceGroupingAnalysisResult::ContainsGantryTilt() -{ - return m_GantryTilt; -} - -void DicomSeriesReader::SliceGroupingAnalysisResult::AddFileToSortedBlock(const std::string& filename) -{ - m_GroupedFiles.push_back( filename ); -} - -void DicomSeriesReader::SliceGroupingAnalysisResult::AddFileToUnsortedBlock(const std::string& filename) -{ - m_UnsortedFiles.push_back( filename ); -} - -void DicomSeriesReader::SliceGroupingAnalysisResult::AddFilesToUnsortedBlock(const StringContainer& filenames) -{ - m_UnsortedFiles.insert( m_UnsortedFiles.end(), filenames.begin(), filenames.end() ); -} - -void DicomSeriesReader::SliceGroupingAnalysisResult::FlagGantryTilt() -{ - m_GantryTilt = true; -} - -void DicomSeriesReader::SliceGroupingAnalysisResult::UndoPrematureGrouping() -{ - assert( !m_GroupedFiles.empty() ); - m_UnsortedFiles.insert( m_UnsortedFiles.begin(), m_GroupedFiles.back() ); - m_GroupedFiles.pop_back(); -} - - - - - - - const DicomSeriesReader::TagToPropertyMapType& DicomSeriesReader::GetDICOMTagsToMITKPropertyMap() { static bool initialized = false; static TagToPropertyMapType dictionary; if (!initialized) { /* Selection criteria: - no sequences because we cannot represent that - nothing animal related (specied, breed registration number), MITK focusses on human medical image processing. - only general attributes so far When extending this, we should make use of a real dictionary (GDCM/DCMTK and lookup the tag names there) */ // Patient module dictionary["0010|0010"] = "dicom.patient.PatientsName"; dictionary["0010|0020"] = "dicom.patient.PatientID"; dictionary["0010|0030"] = "dicom.patient.PatientsBirthDate"; dictionary["0010|0040"] = "dicom.patient.PatientsSex"; dictionary["0010|0032"] = "dicom.patient.PatientsBirthTime"; dictionary["0010|1000"] = "dicom.patient.OtherPatientIDs"; dictionary["0010|1001"] = "dicom.patient.OtherPatientNames"; dictionary["0010|2160"] = "dicom.patient.EthnicGroup"; dictionary["0010|4000"] = "dicom.patient.PatientComments"; dictionary["0012|0062"] = "dicom.patient.PatientIdentityRemoved"; dictionary["0012|0063"] = "dicom.patient.DeIdentificationMethod"; // General Study module dictionary["0020|000d"] = "dicom.study.StudyInstanceUID"; dictionary["0008|0020"] = "dicom.study.StudyDate"; dictionary["0008|0030"] = "dicom.study.StudyTime"; dictionary["0008|0090"] = "dicom.study.ReferringPhysiciansName"; dictionary["0020|0010"] = "dicom.study.StudyID"; dictionary["0008|0050"] = "dicom.study.AccessionNumber"; dictionary["0008|1030"] = "dicom.study.StudyDescription"; dictionary["0008|1048"] = "dicom.study.PhysiciansOfRecord"; dictionary["0008|1060"] = "dicom.study.NameOfPhysicianReadingStudy"; // General Series module dictionary["0008|0060"] = "dicom.series.Modality"; dictionary["0020|000e"] = "dicom.series.SeriesInstanceUID"; dictionary["0020|0011"] = "dicom.series.SeriesNumber"; dictionary["0020|0060"] = "dicom.series.Laterality"; dictionary["0008|0021"] = "dicom.series.SeriesDate"; dictionary["0008|0031"] = "dicom.series.SeriesTime"; dictionary["0008|1050"] = "dicom.series.PerformingPhysiciansName"; dictionary["0018|1030"] = "dicom.series.ProtocolName"; dictionary["0008|103e"] = "dicom.series.SeriesDescription"; dictionary["0008|1070"] = "dicom.series.OperatorsName"; dictionary["0018|0015"] = "dicom.series.BodyPartExamined"; dictionary["0018|5100"] = "dicom.series.PatientPosition"; dictionary["0028|0108"] = "dicom.series.SmallestPixelValueInSeries"; dictionary["0028|0109"] = "dicom.series.LargestPixelValueInSeries"; // VOI LUT module dictionary["0028|1050"] = "dicom.voilut.WindowCenter"; dictionary["0028|1051"] = "dicom.voilut.WindowWidth"; dictionary["0028|1055"] = "dicom.voilut.WindowCenterAndWidthExplanation"; // Image Pixel module dictionary["0028|0004"] = "dicom.pixel.PhotometricInterpretation"; dictionary["0028|0010"] = "dicom.pixel.Rows"; dictionary["0028|0011"] = "dicom.pixel.Columns"; // Image Plane module dictionary["0028|0030"] = "dicom.PixelSpacing"; dictionary["0018|1164"] = "dicom.ImagerPixelSpacing"; initialized = true; } return dictionary; } DataNode::Pointer DicomSeriesReader::LoadDicomSeries(const StringContainer &filenames, bool sort, bool check_4d, bool correctTilt, UpdateCallBackMethod callback, Image::Pointer preLoadedImageBlock) { DataNode::Pointer node = DataNode::New(); if (DicomSeriesReader::LoadDicomSeries(filenames, *node, sort, check_4d, correctTilt, callback, preLoadedImageBlock)) { if( filenames.empty() ) { return NULL; } return node; } else { return NULL; } } bool DicomSeriesReader::LoadDicomSeries( const StringContainer &filenames, DataNode &node, bool sort, bool check_4d, bool correctTilt, UpdateCallBackMethod callback, Image::Pointer preLoadedImageBlock) { if( filenames.empty() ) { MITK_DEBUG << "Calling LoadDicomSeries with empty filename string container. Probably invalid application logic."; node.SetData(NULL); return true; // this is not actually an error but the result is very simple } DcmIoType::Pointer io = DcmIoType::New(); try { if (io->CanReadFile(filenames.front().c_str())) { io->SetFileName(filenames.front().c_str()); io->ReadImageInformation(); - if (io->GetPixelType() == itk::ImageIOBase::SCALAR) + if (io->GetPixelType() == itk::ImageIOBase::SCALAR || + io->GetPixelType() == itk::ImageIOBase::RGB) { - switch (io->GetComponentType()) - { - case DcmIoType::UCHAR: - DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); - break; - case DcmIoType::CHAR: - DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); - break; - case DcmIoType::USHORT: - DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); - break; - case DcmIoType::SHORT: - DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); - break; - case DcmIoType::UINT: - DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); - break; - case DcmIoType::INT: - DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); - break; - case DcmIoType::ULONG: - DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); - break; - case DcmIoType::LONG: - DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); - break; - case DcmIoType::FLOAT: - DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); - break; - case DcmIoType::DOUBLE: - DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); - break; - default: - MITK_ERROR << "Found unsupported DICOM scalar pixel type: (enum value) " << io->GetComponentType(); - } - } - else if (io->GetPixelType() == itk::ImageIOBase::RGB) - { - switch (io->GetComponentType()) - { - case DcmIoType::UCHAR: - DicomSeriesReader::LoadDicom< itk::RGBPixel >(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); - break; - case DcmIoType::CHAR: - DicomSeriesReader::LoadDicom >(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); - break; - case DcmIoType::USHORT: - DicomSeriesReader::LoadDicom >(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); - break; - case DcmIoType::SHORT: - DicomSeriesReader::LoadDicom >(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); - break; - case DcmIoType::UINT: - DicomSeriesReader::LoadDicom >(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); - break; - case DcmIoType::INT: - DicomSeriesReader::LoadDicom >(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); - break; - case DcmIoType::ULONG: - DicomSeriesReader::LoadDicom >(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); - break; - case DcmIoType::LONG: - DicomSeriesReader::LoadDicom >(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); - break; - case DcmIoType::FLOAT: - DicomSeriesReader::LoadDicom >(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); - break; - case DcmIoType::DOUBLE: - DicomSeriesReader::LoadDicom >(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); - break; - default: - MITK_ERROR << "Found unsupported DICOM scalar pixel type: (enum value) " << io->GetComponentType(); - } + LoadDicom(filenames, node, sort, check_4d, correctTilt, callback, preLoadedImageBlock); } if (node.GetData()) { return true; } } } catch(itk::MemoryAllocationError& e) { MITK_ERROR << "Out of memory. Cannot load DICOM series: " << e.what(); } catch(std::exception& e) { MITK_ERROR << "Error encountered when loading DICOM series:" << e.what(); } catch(...) { MITK_ERROR << "Unspecified error encountered when loading DICOM series."; } return false; } bool DicomSeriesReader::IsDicom(const std::string &filename) { DcmIoType::Pointer io = DcmIoType::New(); return io->CanReadFile(filename.c_str()); } bool DicomSeriesReader::IsPhilips3DDicom(const std::string &filename) { DcmIoType::Pointer io = DcmIoType::New(); if (io->CanReadFile(filename.c_str())) { //Look at header Tag 3001,0010 if it is "Philips3D" gdcm::Reader reader; reader.SetFileName(filename.c_str()); reader.Read(); gdcm::DataSet &data_set = reader.GetFile().GetDataSet(); gdcm::StringFilter sf; sf.SetFile(reader.GetFile()); if (data_set.FindDataElement(gdcm::Tag(0x3001, 0x0010)) && (sf.ToString(gdcm::Tag(0x3001, 0x0010)) == "Philips3D ")) { return true; } } return false; } bool DicomSeriesReader::ReadPhilips3DDicom(const std::string &filename, mitk::Image::Pointer output_image) { // Now get PhilipsSpecific Tags gdcm::PixmapReader reader; reader.SetFileName(filename.c_str()); reader.Read(); gdcm::DataSet &data_set = reader.GetFile().GetDataSet(); gdcm::StringFilter sf; sf.SetFile(reader.GetFile()); gdcm::Attribute<0x0028,0x0011> dimTagX; // coloumns || sagittal gdcm::Attribute<0x3001,0x1001, gdcm::VR::UL, gdcm::VM::VM1> dimTagZ; //I have no idea what is VM1. // (Philips specific) // axial gdcm::Attribute<0x0028,0x0010> dimTagY; // rows || coronal gdcm::Attribute<0x0028,0x0008> dimTagT; // how many frames gdcm::Attribute<0x0018,0x602c> spaceTagX; // Spacing in X , unit is "physicalTagx" (usually centimeter) gdcm::Attribute<0x0018,0x602e> spaceTagY; gdcm::Attribute<0x3001,0x1003, gdcm::VR::FD, gdcm::VM::VM1> spaceTagZ; // (Philips specific) gdcm::Attribute<0x0018,0x6024> physicalTagX; // if 3, then spacing params are centimeter gdcm::Attribute<0x0018,0x6026> physicalTagY; gdcm::Attribute<0x3001,0x1002, gdcm::VR::US, gdcm::VM::VM1> physicalTagZ; // (Philips specific) dimTagX.Set(data_set); dimTagY.Set(data_set); dimTagZ.Set(data_set); dimTagT.Set(data_set); spaceTagX.Set(data_set); spaceTagY.Set(data_set); spaceTagZ.Set(data_set); physicalTagX.Set(data_set); physicalTagY.Set(data_set); physicalTagZ.Set(data_set); unsigned int dimX = dimTagX.GetValue(), dimY = dimTagY.GetValue(), dimZ = dimTagZ.GetValue(), dimT = dimTagT.GetValue(), physicalX = physicalTagX.GetValue(), physicalY = physicalTagY.GetValue(), physicalZ = physicalTagZ.GetValue(); float spaceX = spaceTagX.GetValue(), spaceY = spaceTagY.GetValue(), spaceZ = spaceTagZ.GetValue(); if (physicalX == 3) // spacing parameter in cm, have to convert it to mm. spaceX = spaceX * 10; if (physicalY == 3) // spacing parameter in cm, have to convert it to mm. spaceY = spaceY * 10; if (physicalZ == 3) // spacing parameter in cm, have to convert it to mm. spaceZ = spaceZ * 10; // Ok, got all necessary Tags! // Now read Pixeldata (7fe0,0010) X x Y x Z x T Elements const gdcm::Pixmap &pixels = reader.GetPixmap(); gdcm::RAWCodec codec; codec.SetPhotometricInterpretation(gdcm::PhotometricInterpretation::MONOCHROME2); codec.SetPixelFormat(pixels.GetPixelFormat()); codec.SetPlanarConfiguration(0); gdcm::DataElement out; codec.Decode(data_set.GetDataElement(gdcm::Tag(0x7fe0, 0x0010)), out); const gdcm::ByteValue *bv = out.GetByteValue(); const char *new_pixels = bv->GetPointer(); // Create MITK Image + Geometry typedef itk::Image ImageType; //Pixeltype might be different sometimes? Maybe read it out from header ImageType::RegionType myRegion; ImageType::SizeType mySize; ImageType::IndexType myIndex; ImageType::SpacingType mySpacing; ImageType::Pointer imageItk = ImageType::New(); mySpacing[0] = spaceX; mySpacing[1] = spaceY; mySpacing[2] = spaceZ; mySpacing[3] = 1; myIndex[0] = 0; myIndex[1] = 0; myIndex[2] = 0; myIndex[3] = 0; mySize[0] = dimX; mySize[1] = dimY; mySize[2] = dimZ; mySize[3] = dimT; myRegion.SetSize( mySize); myRegion.SetIndex( myIndex ); imageItk->SetSpacing(mySpacing); imageItk->SetRegions( myRegion); imageItk->Allocate(); imageItk->FillBuffer(0); itk::ImageRegionIterator iterator(imageItk, imageItk->GetLargestPossibleRegion()); iterator.GoToBegin(); unsigned long pixCount = 0; unsigned long planeSize = dimX*dimY; unsigned long planeCount = 0; unsigned long timeCount = 0; unsigned long numberOfSlices = dimZ; while (!iterator.IsAtEnd()) { unsigned long adressedPixel = pixCount + (numberOfSlices-1-planeCount)*planeSize // add offset to adress the first pixel of current plane + timeCount*numberOfSlices*planeSize; // add time offset iterator.Set( new_pixels[ adressedPixel ] ); pixCount++; ++iterator; if (pixCount == planeSize) { pixCount = 0; planeCount++; } if (planeCount == numberOfSlices) { planeCount = 0; timeCount++; } if (timeCount == dimT) { break; } } mitk::CastToMitkImage(imageItk, output_image); return true; // actually never returns false yet.. but exception possible } -DicomSeriesReader::GantryTiltInformation::GantryTiltInformation() -: m_ShiftUp(0.0) -, m_ShiftRight(0.0) -, m_ShiftNormal(0.0) -, m_ITKAssumedSliceSpacing(0.0) -, m_NumberOfSlicesApart(1) -{ -} - - -#define doublepoint(x) \ - Point3Dd x; \ - x[0] = x ## f[0]; \ - x[1] = x ## f[1]; \ - x[2] = x ## f[2]; - - -#define doublevector(x) \ - Vector3Dd x; \ - x[0] = x ## f[0]; \ - x[1] = x ## f[1]; \ - x[2] = x ## f[2]; - -DicomSeriesReader::GantryTiltInformation::GantryTiltInformation( - const Point3D& origin1f, const Point3D& origin2f, - const Vector3D& rightf, const Vector3D& upf, - unsigned int numberOfSlicesApart) -: m_ShiftUp(0.0) -, m_ShiftRight(0.0) -, m_ShiftNormal(0.0) -, m_NumberOfSlicesApart(numberOfSlicesApart) -{ - assert(numberOfSlicesApart); - - doublepoint(origin1); - doublepoint(origin2); - doublevector(right); - doublevector(up); - - // determine if slice 1 (imagePosition1 and imageOrientation1) and slice 2 can be in one orthogonal slice stack: - // calculate a line from origin 1, directed along the normal of slice (calculated as the cross product of orientation 1) - // check if this line passes through origin 2 - - /* - Determine if line (imagePosition2 + l * normal) contains imagePosition1. - Done by calculating the distance of imagePosition1 from line (imagePosition2 + l *normal) - - E.g. http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html - - squared distance = | (pointAlongNormal - origin2) x (origin2 - origin1) | ^ 2 - / - |pointAlongNormal - origin2| ^ 2 - - ( x meaning the cross product ) - */ - - Vector3Dd normal = itk::CrossProduct(right, up); - Point3Dd pointAlongNormal = origin2 + normal; - - double numerator = itk::CrossProduct( pointAlongNormal - origin2 , origin2 - origin1 ).GetSquaredNorm(); - double denominator = (pointAlongNormal - origin2).GetSquaredNorm(); - - double distance = sqrt(numerator / denominator); - - if ( distance > 0.001 ) // mitk::eps is too small; 1/1000 of a mm should be enough to detect tilt - { - MITK_DEBUG << " Series seems to contain a tilted (or sheared) geometry"; - MITK_DEBUG << " Distance of expected slice origin from actual slice origin: " << distance; - MITK_DEBUG << " ==> storing this shift for later analysis:"; - MITK_DEBUG << " v right: " << right; - MITK_DEBUG << " v up: " << up; - MITK_DEBUG << " v normal: " << normal; - - Point3Dd projectionRight = projectPointOnLine( origin1, origin2, right ); - Point3Dd projectionNormal = projectPointOnLine( origin1, origin2, normal ); - - m_ShiftRight = (projectionRight - origin2).GetNorm(); - m_ShiftNormal = (projectionNormal - origin2).GetNorm(); - - /* - now also check to which side the image is shifted. - - Calculation e.g. from - http://mathworld.wolfram.com/Point-PlaneDistance.html - */ - - Point3Dd testPoint = origin1; - Vector3Dd planeNormal = up; - - double signedDistance = ( - planeNormal[0] * testPoint[0] - + planeNormal[1] * testPoint[1] - + planeNormal[2] * testPoint[2] - - ( - planeNormal[0] * origin2[0] - + planeNormal[1] * origin2[1] - + planeNormal[2] * origin2[2] - ) - ) - / - sqrt( planeNormal[0] * planeNormal[0] - + planeNormal[1] * planeNormal[1] - + planeNormal[2] * planeNormal[2] - ); - - m_ShiftUp = signedDistance; - - m_ITKAssumedSliceSpacing = (origin2 - origin1).GetNorm(); - // How do we now this is assumed? See header documentation for ITK code references - //double itkAssumedSliceSpacing = sqrt( m_ShiftUp * m_ShiftUp + m_ShiftNormal * m_ShiftNormal ); - - MITK_DEBUG << " shift normal: " << m_ShiftNormal; - MITK_DEBUG << " shift normal assumed by ITK: " << m_ITKAssumedSliceSpacing; - MITK_DEBUG << " shift up: " << m_ShiftUp; - MITK_DEBUG << " shift right: " << m_ShiftRight; - - MITK_DEBUG << " tilt angle (deg): " << atan( m_ShiftUp / m_ShiftNormal ) * 180.0 / 3.1415926535; - } -} - -Point3D -DicomSeriesReader::GantryTiltInformation::projectPointOnLine( Point3Dd p, Point3Dd lineOrigin, Vector3Dd lineDirection ) -{ - /** - See illustration at http://mo.mathematik.uni-stuttgart.de/inhalt/aussage/aussage472/ - - vector(lineOrigin,p) = normal * ( innerproduct((p - lineOrigin),normal) / squared-length(normal) ) - */ - - Vector3Dd lineOriginToP = p - lineOrigin; - double innerProduct = lineOriginToP * lineDirection; - - double factor = innerProduct / lineDirection.GetSquaredNorm(); - Point3Dd projection = lineOrigin + factor * lineDirection; - - return projection; -} - -double -DicomSeriesReader::GantryTiltInformation::GetTiltCorrectedAdditionalSize() const -{ - return fabs(m_ShiftUp); -} - -double -DicomSeriesReader::GantryTiltInformation::GetTiltAngleInDegrees() const -{ - return atan( fabs(m_ShiftUp) / m_ShiftNormal ) * 180.0 / 3.1415926535; -} - -double -DicomSeriesReader::GantryTiltInformation::GetMatrixCoefficientForCorrectionInWorldCoordinates() const -{ - // so many mm need to be shifted per slice! - return m_ShiftUp / static_cast(m_NumberOfSlicesApart); -} - -double -DicomSeriesReader::GantryTiltInformation::GetRealZSpacing() const -{ - return m_ShiftNormal / static_cast(m_NumberOfSlicesApart); -} - - -bool -DicomSeriesReader::GantryTiltInformation::IsSheared() const -{ - return ( fabs(m_ShiftRight) > 0.001 - || fabs(m_ShiftUp) > 0.001); -} - - -bool -DicomSeriesReader::GantryTiltInformation::IsRegularGantryTilt() const -{ - return ( fabs(m_ShiftRight) < 0.001 - && fabs(m_ShiftUp) > 0.001); -} - - std::string DicomSeriesReader::ConstCharStarToString(const char* s) { return s ? std::string(s) : std::string(); } bool DicomSeriesReader::DICOMStringToSpacing(const std::string& s, float& spacingX, float& spacingY) { bool successful = false; std::istringstream spacingReader(s); std::string spacing; if ( std::getline( spacingReader, spacing, '\\' ) ) { spacingY = atof( spacing.c_str() ); if ( std::getline( spacingReader, spacing, '\\' ) ) { spacingX = atof( spacing.c_str() ); successful = true; } } return successful; } Point3D DicomSeriesReader::DICOMStringToPoint3D(const std::string& s, bool& successful) { Point3D p; successful = true; std::istringstream originReader(s); std::string coordinate; unsigned int dim(0); while( std::getline( originReader, coordinate, '\\' ) && dim < 3) { p[dim++]= atof(coordinate.c_str()); } if (dim && dim != 3) { successful = false; MITK_ERROR << "Reader implementation made wrong assumption on tag (0020,0032). Found " << dim << " instead of 3 values."; } else if (dim == 0) { successful = false; p.Fill(0.0); // assume default (0,0,0) } return p; } void DicomSeriesReader::DICOMStringToOrientationVectors(const std::string& s, Vector3D& right, Vector3D& up, bool& successful) { successful = true; std::istringstream orientationReader(s); std::string coordinate; unsigned int dim(0); while( std::getline( orientationReader, coordinate, '\\' ) && dim < 6 ) { if (dim<3) { right[dim++] = atof(coordinate.c_str()); } else { up[dim++ - 3] = atof(coordinate.c_str()); } } if (dim && dim != 6) { successful = false; MITK_ERROR << "Reader implementation made wrong assumption on tag (0020,0037). Found " << dim << " instead of 6 values."; } else if (dim == 0) { // fill with defaults right.Fill(0.0); right[0] = 1.0; up.Fill(0.0); up[1] = 1.0; successful = false; } } DicomSeriesReader::SliceGroupingAnalysisResult DicomSeriesReader::AnalyzeFileForITKImageSeriesReaderSpacingAssumption( const StringContainer& files, bool groupImagesWithGantryTilt, const gdcm::Scanner::MappingType& tagValueMappings_) { // result.first = files that fit ITK's assumption // result.second = files that do not fit, should be run through AnalyzeFileForITKImageSeriesReaderSpacingAssumption() again SliceGroupingAnalysisResult result; // we const_cast here, because I could not use a map.at(), which would make the code much more readable gdcm::Scanner::MappingType& tagValueMappings = const_cast(tagValueMappings_); const gdcm::Tag tagImagePositionPatient(0x0020,0x0032); // Image Position (Patient) const gdcm::Tag tagImageOrientation(0x0020, 0x0037); // Image Orientation const gdcm::Tag tagGantryTilt(0x0018, 0x1120); // gantry tilt Vector3D fromFirstToSecondOrigin; fromFirstToSecondOrigin.Fill(0.0); bool fromFirstToSecondOriginInitialized(false); Point3D thisOrigin; thisOrigin.Fill(0.0f); Point3D lastOrigin; lastOrigin.Fill(0.0f); Point3D lastDifferentOrigin; lastDifferentOrigin.Fill(0.0f); bool lastOriginInitialized(false); MITK_DEBUG << "--------------------------------------------------------------------------------"; MITK_DEBUG << "Analyzing files for z-spacing assumption of ITK's ImageSeriesReader (group tilted: " << groupImagesWithGantryTilt << ")"; unsigned int fileIndex(0); for (StringContainer::const_iterator fileIter = files.begin(); fileIter != files.end(); ++fileIter, ++fileIndex) { bool fileFitsIntoPattern(false); std::string thisOriginString; // Read tag value into point3D. PLEASE replace this by appropriate GDCM code if you figure out how to do that thisOriginString = ConstCharStarToString( tagValueMappings[fileIter->c_str()][tagImagePositionPatient] ); if (thisOriginString.empty()) { // don't let such files be in a common group. Everything without position information will be loaded as a single slice: // with standard DICOM files this can happen to: CR, DX, SC MITK_DEBUG << " ==> Sort away " << *fileIter << " for later analysis (no position information)"; // we already have one occupying this position if ( result.GetBlockFilenames().empty() ) // nothing WITH position information yet { // ==> this is a group of its own, stop processing, come back later result.AddFileToSortedBlock( *fileIter ); StringContainer remainingFiles; remainingFiles.insert( remainingFiles.end(), fileIter+1, files.end() ); result.AddFilesToUnsortedBlock( remainingFiles ); fileFitsIntoPattern = false; break; // no files anymore } else { // ==> this does not match, consider later result.AddFileToUnsortedBlock( *fileIter ); fileFitsIntoPattern = false; continue; // next file } } bool ignoredConversionError(-42); // hard to get here, no graceful way to react thisOrigin = DICOMStringToPoint3D( thisOriginString, ignoredConversionError ); MITK_DEBUG << " " << fileIndex << " " << *fileIter << " at " /* << thisOriginString */ << "(" << thisOrigin[0] << "," << thisOrigin[1] << "," << thisOrigin[2] << ")"; if ( lastOriginInitialized && (thisOrigin == lastOrigin) ) { MITK_DEBUG << " ==> Sort away " << *fileIter << " for separate time step"; // we already have one occupying this position result.AddFileToUnsortedBlock( *fileIter ); fileFitsIntoPattern = false; } else { if (!fromFirstToSecondOriginInitialized && lastOriginInitialized) // calculate vector as soon as possible when we get a new position { fromFirstToSecondOrigin = thisOrigin - lastDifferentOrigin; fromFirstToSecondOriginInitialized = true; // Here we calculate if this slice and the previous one are well aligned, // i.e. we test if the previous origin is on a line through the current // origin, directed into the normal direction of the current slice. // If this is NOT the case, then we have a data set with a TILTED GANTRY geometry, // which cannot be simply loaded into a single mitk::Image at the moment. // For this case, we flag this finding in the result and DicomSeriesReader // can correct for that later. Vector3D right; right.Fill(0.0); Vector3D up; right.Fill(0.0); // might be down as well, but it is just a name at this point DICOMStringToOrientationVectors( tagValueMappings[fileIter->c_str()][tagImageOrientation], right, up, ignoredConversionError ); GantryTiltInformation tiltInfo( lastDifferentOrigin, thisOrigin, right, up, 1 ); if ( tiltInfo.IsSheared() ) // mitk::eps is too small; 1/1000 of a mm should be enough to detect tilt { /* optimistic approach, accepting gantry tilt: save file for later, check all further files */ // at this point we have TWO slices analyzed! if they are the only two files, we still split, because there is no third to verify our tilting assumption. // later with a third being available, we must check if the initial tilting vector is still valid. if yes, continue. // if NO, we need to split the already sorted part (result.first) and the currently analyzed file (*fileIter) // tell apart gantry tilt from overall skewedness // sort out irregularly sheared slices, that IS NOT tilting if ( groupImagesWithGantryTilt && tiltInfo.IsRegularGantryTilt() ) { // check if this is at least roughly the same angle as recorded in DICOM tags if ( tagValueMappings[fileIter->c_str()].find(tagGantryTilt) != tagValueMappings[fileIter->c_str()].end() ) { // read value, compare to calculated angle std::string tiltStr = ConstCharStarToString( tagValueMappings[fileIter->c_str()][tagGantryTilt] ); double angle = atof(tiltStr.c_str()); MITK_DEBUG << "Comparing recorded tilt angle " << angle << " against calculated value " << tiltInfo.GetTiltAngleInDegrees(); // TODO we probably want the signs correct, too (that depends: this is just a rough check, nothing serious) if ( fabs(angle) - tiltInfo.GetTiltAngleInDegrees() > 0.25) { result.AddFileToUnsortedBlock( *fileIter ); // sort away for further analysis fileFitsIntoPattern = false; } else // tilt angle from header is less than 0.25 degrees different from what we calculated, assume this is fine { result.FlagGantryTilt(); result.AddFileToSortedBlock(*fileIter); // this file is good for current block fileFitsIntoPattern = true; } } else // we cannot check the calculated tilt angle against the one from the dicom header (so we assume we are right) { result.FlagGantryTilt(); result.AddFileToSortedBlock(*fileIter); // this file is good for current block fileFitsIntoPattern = true; } } else // caller does not want tilt compensation OR shearing is more complicated than tilt { result.AddFileToUnsortedBlock( *fileIter ); // sort away for further analysis fileFitsIntoPattern = false; } } else // not sheared { result.AddFileToSortedBlock(*fileIter); // this file is good for current block fileFitsIntoPattern = true; } } else if (fromFirstToSecondOriginInitialized) // we already know the offset between slices { Point3D assumedOrigin = lastDifferentOrigin + fromFirstToSecondOrigin; Vector3D originError = assumedOrigin - thisOrigin; double norm = originError.GetNorm(); double toleratedError(0.005); // max. 1/10mm error when measurement crosses 20 slices in z direction if (norm > toleratedError) { MITK_DEBUG << " File does not fit into the inter-slice distance pattern (diff = " << norm << ", allowed " << toleratedError << ")."; MITK_DEBUG << " Expected position (" << assumedOrigin[0] << "," << assumedOrigin[1] << "," << assumedOrigin[2] << "), got position (" << thisOrigin[0] << "," << thisOrigin[1] << "," << thisOrigin[2] << ")"; MITK_DEBUG << " ==> Sort away " << *fileIter << " for later analysis"; // At this point we know we deviated from the expectation of ITK's ImageSeriesReader // We split the input file list at this point, i.e. all files up to this one (excluding it) // are returned as group 1, the remaining files (including the faulty one) are group 2 /* Optimistic approach: check if any of the remaining slices fits in */ result.AddFileToUnsortedBlock( *fileIter ); // sort away for further analysis fileFitsIntoPattern = false; } else { result.AddFileToSortedBlock(*fileIter); // this file is good for current block fileFitsIntoPattern = true; } } else // this should be the very first slice { result.AddFileToSortedBlock(*fileIter); // this file is good for current block fileFitsIntoPattern = true; } } // record current origin for reference in later iterations if ( !lastOriginInitialized || ( fileFitsIntoPattern && (thisOrigin != lastOrigin) ) ) { lastDifferentOrigin = thisOrigin; } lastOrigin = thisOrigin; lastOriginInitialized = true; } if ( result.ContainsGantryTilt() ) { // check here how many files were grouped. // IF it was only two files AND we assume tiltedness (e.g. save "distance") // THEN we would want to also split the two previous files (simple) because // we don't have any reason to assume they belong together if ( result.GetBlockFilenames().size() == 2 ) { result.UndoPrematureGrouping(); } } return result; } DicomSeriesReader::FileNamesGrouping DicomSeriesReader::GetSeries(const StringContainer& files, bool groupImagesWithGantryTilt, const StringContainer &restrictions) { return GetSeries(files, true, groupImagesWithGantryTilt, restrictions); } DicomSeriesReader::FileNamesGrouping DicomSeriesReader::GetSeries(const StringContainer& files, bool sortTo3DPlust, bool groupImagesWithGantryTilt, const StringContainer& /*restrictions*/) { /** assumption about this method: returns a map of uid-like-key --> list(filename) each entry should contain filenames that have images of same - series instance uid (automatically done by GDCMSeriesFileNames - 0020,0037 image orientation (patient) - 0028,0030 pixel spacing (x,y) - 0018,0050 slice thickness */ // use GDCM directly, itk::GDCMSeriesFileNames does not work with GDCM 2 // PART I: scan files for sorting relevant DICOM tags, // separate images that differ in any of those // attributes (they cannot possibly form a 3D block) // scan for relevant tags in dicom files gdcm::Scanner scanner; const gdcm::Tag tagSOPClassUID(0x0008, 0x0016); // SOP class UID scanner.AddTag( tagSOPClassUID ); const gdcm::Tag tagSeriesInstanceUID(0x0020,0x000e); // Series Instance UID scanner.AddTag( tagSeriesInstanceUID ); const gdcm::Tag tagImageOrientation(0x0020, 0x0037); // image orientation scanner.AddTag( tagImageOrientation ); const gdcm::Tag tagPixelSpacing(0x0028, 0x0030); // pixel spacing scanner.AddTag( tagPixelSpacing ); const gdcm::Tag tagImagerPixelSpacing(0x0018, 0x1164); // imager pixel spacing scanner.AddTag( tagImagerPixelSpacing ); const gdcm::Tag tagSliceThickness(0x0018, 0x0050); // slice thickness scanner.AddTag( tagSliceThickness ); const gdcm::Tag tagNumberOfRows(0x0028, 0x0010); // number rows scanner.AddTag( tagNumberOfRows ); const gdcm::Tag tagNumberOfColumns(0x0028, 0x0011); // number cols scanner.AddTag( tagNumberOfColumns ); const gdcm::Tag tagGantryTilt(0x0018, 0x1120); // gantry tilt scanner.AddTag( tagGantryTilt ); const gdcm::Tag tagModality(0x0008, 0x0060); // modality scanner.AddTag( tagModality ); const gdcm::Tag tagNumberOfFrames(0x0028, 0x0008); // number of frames scanner.AddTag( tagNumberOfFrames ); // additional tags read in this scan to allow later analysis // THESE tag are not used for initial separating of files const gdcm::Tag tagImagePositionPatient(0x0020,0x0032); // Image Position (Patient) scanner.AddTag( tagImagePositionPatient ); // TODO add further restrictions from arguments (when anybody asks for it) FileNamesGrouping result; // let GDCM scan files if ( !scanner.Scan( files ) ) { MITK_ERROR << "gdcm::Scanner failed when scanning " << files.size() << " input files."; return result; } // assign files IDs that will separate them for loading into image blocks for (gdcm::Scanner::ConstIterator fileIter = scanner.Begin(); fileIter != scanner.End(); ++fileIter) { if ( std::string(fileIter->first).empty() ) continue; // TODO understand why Scanner has empty string entries if ( std::string(fileIter->first) == std::string("DICOMDIR") ) continue; /* sort out multi-frame if ( scanner.GetValue( fileIter->first , tagNumberOfFrames ) ) { MITK_INFO << "Ignoring " << fileIter->first << " because we cannot handle multi-frame images."; continue; } */ // we const_cast here, because I could not use a map.at() function in CreateMoreUniqueSeriesIdentifier. // doing the same thing with find would make the code less readable. Since we forget the Scanner results // anyway after this function, we can simply tolerate empty map entries introduced by bad operator[] access std::string moreUniqueSeriesId = CreateMoreUniqueSeriesIdentifier( const_cast(fileIter->second) ); result[ moreUniqueSeriesId ].AddFile( fileIter->first ); } // PART II: sort slices spatially (or at least consistently if this is NOT possible, see method) for ( FileNamesGrouping::const_iterator groupIter = result.begin(); groupIter != result.end(); ++groupIter ) { try { result[ groupIter->first ] = ImageBlockDescriptor( SortSeriesSlices( groupIter->second.GetFilenames() ) ); // sort each slice group spatially } catch(...) { MITK_ERROR << "Caught something."; } } // PART III: analyze pre-sorted images for valid blocks (i.e. blocks of equal z-spacing), // separate into multiple blocks if necessary. // // Analysis performs the following steps: // * imitate itk::ImageSeriesReader: use the distance between the first two images as z-spacing // * check what images actually fulfill ITK's z-spacing assumption // * separate all images that fail the test into new blocks, re-iterate analysis for these blocks // * this includes images which DO NOT PROVIDE spatial information, i.e. all images w/o ImagePositionPatient will be loaded separately FileNamesGrouping groupsOf3DPlusTBlocks; // final result of this function for ( FileNamesGrouping::const_iterator groupIter = result.begin(); groupIter != result.end(); ++groupIter ) { FileNamesGrouping groupsOf3DBlocks; // intermediate result for only this group(!) std::map mapOf3DBlockAnalysisResults; StringContainer filesStillToAnalyze = groupIter->second.GetFilenames(); std::string groupUID = groupIter->first; unsigned int subgroup(0); MITK_DEBUG << "Analyze group " << groupUID; while (!filesStillToAnalyze.empty()) // repeat until all files are grouped somehow { SliceGroupingAnalysisResult analysisResult = AnalyzeFileForITKImageSeriesReaderSpacingAssumption( filesStillToAnalyze, groupImagesWithGantryTilt, scanner.GetMappings() ); // enhance the UID for additional groups std::stringstream newGroupUID; newGroupUID << groupUID << '.' << subgroup; ImageBlockDescriptor thisBlock( analysisResult.GetBlockFilenames() ); std::string firstFileInBlock = thisBlock.GetFilenames().front(); thisBlock.SetImageBlockUID( newGroupUID.str() ); thisBlock.SetSeriesInstanceUID( DicomSeriesReader::ConstCharStarToString( scanner.GetValue( firstFileInBlock.c_str(), tagSeriesInstanceUID ) ) ); thisBlock.SetHasGantryTiltCorrected( analysisResult.ContainsGantryTilt() ); thisBlock.SetSOPClassUID( DicomSeriesReader::ConstCharStarToString( scanner.GetValue( firstFileInBlock.c_str(), tagSOPClassUID ) ) ); thisBlock.SetNumberOfFrames( ConstCharStarToString( scanner.GetValue( firstFileInBlock.c_str(), tagNumberOfFrames ) ) ); thisBlock.SetModality( DicomSeriesReader::ConstCharStarToString( scanner.GetValue( firstFileInBlock.c_str(), tagModality ) ) ); thisBlock.SetPixelSpacingInformation( DicomSeriesReader::ConstCharStarToString( scanner.GetValue( firstFileInBlock.c_str(), tagPixelSpacing ) ), DicomSeriesReader::ConstCharStarToString( scanner.GetValue( firstFileInBlock.c_str(), tagImagerPixelSpacing ) ) ); thisBlock.SetHasMultipleTimePoints( false ); groupsOf3DBlocks[ newGroupUID.str() ] = thisBlock; //MITK_DEBUG << "Result: sorted 3D group " << newGroupUID.str() << " with " << groupsOf3DBlocks[ newGroupUID.str() ].GetFilenames().size() << " files"; MITK_DEBUG << "Result: sorted 3D group with " << groupsOf3DBlocks[ newGroupUID.str() ].GetFilenames().size() << " files"; StringContainer debugOutputFiles = analysisResult.GetBlockFilenames(); for (StringContainer::const_iterator siter = debugOutputFiles.begin(); siter != debugOutputFiles.end(); ++siter) MITK_DEBUG << " IN " << *siter; ++subgroup; filesStillToAnalyze = analysisResult.GetUnsortedFilenames(); // remember what needs further analysis for (StringContainer::const_iterator siter = filesStillToAnalyze.begin(); siter != filesStillToAnalyze.end(); ++siter) MITK_DEBUG << " OUT " << *siter; } // end of grouping, now post-process groups // PART IV: attempt to group blocks to 3D+t blocks if requested // inspect entries of groupsOf3DBlocks // - if number of files is identical to previous entry, collect for 3D+t block // - as soon as number of files changes from previous entry, record collected blocks as 3D+t block, start a new one, continue // decide whether or not to group 3D blocks into 3D+t blocks where possible if ( !sortTo3DPlust ) { // copy 3D blocks to output groupsOf3DPlusTBlocks.insert( groupsOf3DBlocks.begin(), groupsOf3DBlocks.end() ); } else { // sort 3D+t (as described in "PART IV") MITK_DEBUG << "================================================================================"; MITK_DEBUG << "3D+t analysis:"; unsigned int numberOfFilesInPreviousBlock(0); std::string previousBlockKey; for ( FileNamesGrouping::const_iterator block3DIter = groupsOf3DBlocks.begin(); block3DIter != groupsOf3DBlocks.end(); ++block3DIter ) { unsigned int numberOfFilesInThisBlock = block3DIter->second.GetFilenames().size(); std::string thisBlockKey = block3DIter->first; if (numberOfFilesInPreviousBlock == 0) { numberOfFilesInPreviousBlock = numberOfFilesInThisBlock; groupsOf3DPlusTBlocks[thisBlockKey] = block3DIter->second; MITK_DEBUG << " 3D+t group " << thisBlockKey; previousBlockKey = thisBlockKey; } else { bool identicalOrigins; try { // check whether this and the previous block share a comon origin // TODO should be safe, but a little try/catch or other error handling wouldn't hurt const char *origin_value = scanner.GetValue( groupsOf3DBlocks[thisBlockKey].GetFilenames().front().c_str(), tagImagePositionPatient ), *previous_origin_value = scanner.GetValue( groupsOf3DBlocks[previousBlockKey].GetFilenames().front().c_str(), tagImagePositionPatient ), *destination_value = scanner.GetValue( groupsOf3DBlocks[thisBlockKey].GetFilenames().back().c_str(), tagImagePositionPatient ), *previous_destination_value = scanner.GetValue( groupsOf3DBlocks[previousBlockKey].GetFilenames().back().c_str(), tagImagePositionPatient ); if (!origin_value || !previous_origin_value || !destination_value || !previous_destination_value) { identicalOrigins = false; } else { std::string thisOriginString = ConstCharStarToString( origin_value ); std::string previousOriginString = ConstCharStarToString( previous_origin_value ); // also compare last origin, because this might differ if z-spacing is different std::string thisDestinationString = ConstCharStarToString( destination_value ); std::string previousDestinationString = ConstCharStarToString( previous_destination_value ); identicalOrigins = ( (thisOriginString == previousOriginString) && (thisDestinationString == previousDestinationString) ); } } catch(...) { identicalOrigins = false; } if (identicalOrigins && (numberOfFilesInPreviousBlock == numberOfFilesInThisBlock)) { // group with previous block groupsOf3DPlusTBlocks[previousBlockKey].AddFiles( block3DIter->second.GetFilenames() ); groupsOf3DPlusTBlocks[previousBlockKey].SetHasMultipleTimePoints(true); MITK_DEBUG << " --> group enhanced with another timestep"; } else { // start a new block groupsOf3DPlusTBlocks[thisBlockKey] = block3DIter->second; int numberOfTimeSteps = groupsOf3DPlusTBlocks[previousBlockKey].GetFilenames().size() / numberOfFilesInPreviousBlock; MITK_DEBUG << " ==> group closed with " << numberOfTimeSteps << " time steps"; previousBlockKey = thisBlockKey; MITK_DEBUG << " 3D+t group " << thisBlockKey << " started"; } } numberOfFilesInPreviousBlock = numberOfFilesInThisBlock; } } } MITK_DEBUG << "================================================================================"; MITK_DEBUG << "Summary: "; for ( FileNamesGrouping::const_iterator groupIter = groupsOf3DPlusTBlocks.begin(); groupIter != groupsOf3DPlusTBlocks.end(); ++groupIter ) { ImageBlockDescriptor block = groupIter->second; MITK_DEBUG << " " << block.GetFilenames().size() << " '" << block.GetModality() << "' images (" << block.GetSOPClassUIDAsString() << ") in volume " << block.GetImageBlockUID(); MITK_DEBUG << " (gantry tilt : " << (block.HasGantryTiltCorrected()?"Yes":"No") << "; " "pixel spacing : " << PixelSpacingInterpretationToString( block.GetPixelSpacingType() ) << "; " "3D+t: " << (block.HasMultipleTimePoints()?"Yes":"No") << "; " "reader support: " << ReaderImplementationLevelToString( block.GetReaderImplementationLevel() ) << ")"; StringContainer debugOutputFiles = block.GetFilenames(); for (StringContainer::const_iterator siter = debugOutputFiles.begin(); siter != debugOutputFiles.end(); ++siter) MITK_DEBUG << " F " << *siter; } MITK_DEBUG << "================================================================================"; return groupsOf3DPlusTBlocks; } DicomSeriesReader::FileNamesGrouping DicomSeriesReader::GetSeries(const std::string &dir, bool groupImagesWithGantryTilt, const StringContainer &restrictions) { gdcm::Directory directoryLister; directoryLister.Load( dir.c_str(), false ); // non-recursive return GetSeries(directoryLister.GetFilenames(), groupImagesWithGantryTilt, restrictions); } std::string DicomSeriesReader::CreateSeriesIdentifierPart( gdcm::Scanner::TagToValue& tagValueMap, const gdcm::Tag& tag ) { std::string result; try { result = IDifyTagValue( tagValueMap[ tag ] ? tagValueMap[ tag ] : std::string("") ); } catch (std::exception&) { // we are happy with even nothing, this will just group images of a series //MITK_WARN << "Could not access tag " << tag << ": " << e.what(); } return result; } std::string DicomSeriesReader::CreateMoreUniqueSeriesIdentifier( gdcm::Scanner::TagToValue& tagValueMap ) { const gdcm::Tag tagSeriesInstanceUID(0x0020,0x000e); // Series Instance UID const gdcm::Tag tagImageOrientation(0x0020, 0x0037); // image orientation const gdcm::Tag tagPixelSpacing(0x0028, 0x0030); // pixel spacing const gdcm::Tag tagImagerPixelSpacing(0x0018, 0x1164); // imager pixel spacing const gdcm::Tag tagSliceThickness(0x0018, 0x0050); // slice thickness const gdcm::Tag tagNumberOfRows(0x0028, 0x0010); // number rows const gdcm::Tag tagNumberOfColumns(0x0028, 0x0011); // number cols const gdcm::Tag tagNumberOfFrames(0x0028, 0x0008); // number of frames const char* tagSeriesInstanceUid = tagValueMap[tagSeriesInstanceUID]; if (!tagSeriesInstanceUid) { mitkThrow() << "CreateMoreUniqueSeriesIdentifier() could not access series instance UID. Something is seriously wrong with this image, so stopping here."; } std::string constructedID = tagSeriesInstanceUid; constructedID += CreateSeriesIdentifierPart( tagValueMap, tagNumberOfRows ); constructedID += CreateSeriesIdentifierPart( tagValueMap, tagNumberOfColumns ); constructedID += CreateSeriesIdentifierPart( tagValueMap, tagPixelSpacing ); constructedID += CreateSeriesIdentifierPart( tagValueMap, tagImagerPixelSpacing ); constructedID += CreateSeriesIdentifierPart( tagValueMap, tagSliceThickness ); constructedID += CreateSeriesIdentifierPart( tagValueMap, tagNumberOfFrames ); // be a bit tolerant for orienatation, let only the first few digits matter (http://bugs.mitk.org/show_bug.cgi?id=12263) // NOT constructedID += CreateSeriesIdentifierPart( tagValueMap, tagImageOrientation ); if (tagValueMap.find(tagImageOrientation) != tagValueMap.end()) { bool conversionError(false); Vector3D right; right.Fill(0.0); Vector3D up; right.Fill(0.0); DICOMStringToOrientationVectors( tagValueMap[tagImageOrientation], right, up, conversionError ); //string newstring sprintf(simplifiedOrientationString, "%.3f\\%.3f\\%.3f\\%.3f\\%.3f\\%.3f", right[0], right[1], right[2], up[0], up[1], up[2]); std::ostringstream ss; ss.setf(std::ios::fixed, std::ios::floatfield); ss.precision(5); ss << right[0] << "\\" << right[1] << "\\" << right[2] << "\\" << up[0] << "\\" << up[1] << "\\" << up[2]; std::string simplifiedOrientationString(ss.str()); constructedID += IDifyTagValue( simplifiedOrientationString ); } constructedID.resize( constructedID.length() - 1 ); // cut of trailing '.' return constructedID; } std::string DicomSeriesReader::IDifyTagValue(const std::string& value) { std::string IDifiedValue( value ); if (value.empty()) throw std::logic_error("IDifyTagValue() illegaly called with empty tag value"); // Eliminate non-alnum characters, including whitespace... // that may have been introduced by concats. for(std::size_t i=0; i= 'a' && IDifiedValue[i] <= 'z') || (IDifiedValue[i] >= '0' && IDifiedValue[i] <= '9') || (IDifiedValue[i] >= 'A' && IDifiedValue[i] <= 'Z'))) { IDifiedValue.erase(i, 1); } } IDifiedValue += "."; return IDifiedValue; } DicomSeriesReader::StringContainer DicomSeriesReader::GetSeries(const std::string &dir, const std::string &series_uid, bool groupImagesWithGantryTilt, const StringContainer &restrictions) { FileNamesGrouping allSeries = GetSeries(dir, groupImagesWithGantryTilt, restrictions); StringContainer resultingFileList; for ( FileNamesGrouping::const_iterator idIter = allSeries.begin(); idIter != allSeries.end(); ++idIter ) { if ( idIter->first.find( series_uid ) == 0 ) // this ID starts with given series_uid { return idIter->second.GetFilenames(); } } return resultingFileList; } DicomSeriesReader::StringContainer DicomSeriesReader::SortSeriesSlices(const StringContainer &unsortedFilenames) { /* we CAN expect a group of equal - series instance uid - image orientation - pixel spacing - imager pixel spacing - slice thickness - number of rows/columns (each piece of information except the rows/columns might be missing) sorting with GdcmSortFunction tries its best by sorting by spatial position and more hints (acquisition number, acquisition time, trigger time) but will always produce a sorting by falling back to SOP Instance UID. */ gdcm::Sorter sorter; sorter.SetSortFunction(DicomSeriesReader::GdcmSortFunction); try { if (sorter.Sort(unsortedFilenames)) { return sorter.GetFilenames(); } else { MITK_WARN << "Sorting error. Leaving series unsorted."; return unsortedFilenames; } } catch(std::logic_error&) { MITK_WARN << "Sorting error. Leaving series unsorted."; return unsortedFilenames; } } bool DicomSeriesReader::GdcmSortFunction(const gdcm::DataSet &ds1, const gdcm::DataSet &ds2) { // This method MUST accept missing location and position information (and all else, too) // because we cannot rely on anything // (restriction on the sentence before: we have to provide consistent sorting, so we // rely on the minimum information all DICOM files need to provide: SOP Instance UID) /* we CAN expect a group of equal - series instance uid - image orientation - pixel spacing - imager pixel spacing - slice thickness - number of rows/columns */ static const gdcm::Tag tagImagePositionPatient(0x0020,0x0032); // Image Position (Patient) static const gdcm::Tag tagImageOrientation(0x0020, 0x0037); // Image Orientation // see if we have Image Position and Orientation if ( ds1.FindDataElement(tagImagePositionPatient) && ds1.FindDataElement(tagImageOrientation) && ds2.FindDataElement(tagImagePositionPatient) && ds2.FindDataElement(tagImageOrientation) ) { gdcm::Attribute<0x0020,0x0032> image_pos1; // Image Position (Patient) gdcm::Attribute<0x0020,0x0037> image_orientation1; // Image Orientation (Patient) image_pos1.Set(ds1); image_orientation1.Set(ds1); gdcm::Attribute<0x0020,0x0032> image_pos2; gdcm::Attribute<0x0020,0x0037> image_orientation2; image_pos2.Set(ds2); image_orientation2.Set(ds2); /* we tolerate very small differences in image orientation, since we got to know about acquisitions where these values change across a single series (7th decimal digit) (http://bugs.mitk.org/show_bug.cgi?id=12263) still, we want to check if our assumption of 'almost equal' orientations is valid */ for (unsigned int dim = 0; dim < 6; ++dim) { if ( fabs(image_orientation2[dim] - image_orientation1[dim]) > 0.0001 ) { MITK_ERROR << "Dicom images have different orientations."; throw std::logic_error("Dicom images have different orientations. Call GetSeries() first to separate images."); } } double normal[3]; normal[0] = image_orientation1[1] * image_orientation1[5] - image_orientation1[2] * image_orientation1[4]; normal[1] = image_orientation1[2] * image_orientation1[3] - image_orientation1[0] * image_orientation1[5]; normal[2] = image_orientation1[0] * image_orientation1[4] - image_orientation1[1] * image_orientation1[3]; double dist1 = 0.0, dist2 = 0.0; // this computes the distance from world origin (0,0,0) ALONG THE NORMAL of the image planes for (unsigned char i = 0u; i < 3u; ++i) { dist1 += normal[i] * image_pos1[i]; dist2 += normal[i] * image_pos2[i]; } // if we can sort by just comparing the distance, we do exactly that if ( fabs(dist1 - dist2) >= mitk::eps) { // default: compare position return dist1 < dist2; } else // we need to check more properties to distinguish slices { // try to sort by Acquisition Number static const gdcm::Tag tagAcquisitionNumber(0x0020, 0x0012); if (ds1.FindDataElement(tagAcquisitionNumber) && ds2.FindDataElement(tagAcquisitionNumber)) { gdcm::Attribute<0x0020,0x0012> acquisition_number1; // Acquisition number gdcm::Attribute<0x0020,0x0012> acquisition_number2; acquisition_number1.Set(ds1); acquisition_number2.Set(ds2); if (acquisition_number1 != acquisition_number2) { return acquisition_number1 < acquisition_number2; } else // neither position nor acquisition number are good for sorting, so check more { // try to sort by Acquisition Time static const gdcm::Tag tagAcquisitionTime(0x0008, 0x0032); if (ds1.FindDataElement(tagAcquisitionTime) && ds2.FindDataElement(tagAcquisitionTime)) { gdcm::Attribute<0x0008,0x0032> acquisition_time1; // Acquisition time gdcm::Attribute<0x0008,0x0032> acquisition_time2; acquisition_time1.Set(ds1); acquisition_time2.Set(ds2); if (acquisition_time1 != acquisition_time2) { return acquisition_time1 < acquisition_time2; } else // we gave up on image position, acquisition number and acquisition time now { // let's try trigger time static const gdcm::Tag tagTriggerTime(0x0018, 0x1060); if (ds1.FindDataElement(tagTriggerTime) && ds2.FindDataElement(tagTriggerTime)) { gdcm::Attribute<0x0018,0x1060> trigger_time1; // Trigger time gdcm::Attribute<0x0018,0x1060> trigger_time2; trigger_time1.Set(ds1); trigger_time2.Set(ds2); if (trigger_time1 != trigger_time2) { return trigger_time1 < trigger_time2; } // ELSE! // for this and many previous ifs we fall through if nothing lets us sort } // . } // . } // . } } } } // . // LAST RESORT: all valuable information for sorting is missing. // Sort by some meaningless but unique identifiers to satisfy the sort function static const gdcm::Tag tagSOPInstanceUID(0x0008, 0x0018); if (ds1.FindDataElement(tagSOPInstanceUID) && ds2.FindDataElement(tagSOPInstanceUID)) { MITK_DEBUG << "Dicom images are missing attributes for a meaningful sorting, falling back to SOP instance UID comparison."; gdcm::Attribute<0x0008,0x0018> SOPInstanceUID1; // SOP instance UID is mandatory and unique gdcm::Attribute<0x0008,0x0018> SOPInstanceUID2; SOPInstanceUID1.Set(ds1); SOPInstanceUID2.Set(ds2); return SOPInstanceUID1 < SOPInstanceUID2; } else { // no DICOM file should really come down here, this should only be reached with unskillful and unlucky manipulation of files std::string error_message("Malformed DICOM images, which do not even contain a SOP Instance UID."); MITK_ERROR << error_message; throw std::logic_error( error_message ); } } std::string DicomSeriesReader::GetConfigurationString() { std::stringstream configuration; configuration << "MITK_USE_GDCMIO: "; configuration << "true"; configuration << "\n"; configuration << "GDCM_VERSION: "; #ifdef GDCM_MAJOR_VERSION configuration << GDCM_VERSION; #endif //configuration << "\n"; return configuration.str(); } void DicomSeriesReader::CopyMetaDataToImageProperties(StringContainer filenames, const gdcm::Scanner::MappingType &tagValueMappings_, DcmIoType *io, const ImageBlockDescriptor& blockInfo, Image *image) { std::list imageBlock; imageBlock.push_back(filenames); CopyMetaDataToImageProperties(imageBlock, tagValueMappings_, io, blockInfo, image); } void DicomSeriesReader::CopyMetaDataToImageProperties( std::list imageBlock, const gdcm::Scanner::MappingType& tagValueMappings_, DcmIoType* io, const ImageBlockDescriptor& blockInfo, Image* image) { if (!io || !image) return; StringLookupTable filesForSlices; StringLookupTable sliceLocationForSlices; StringLookupTable instanceNumberForSlices; StringLookupTable SOPInstanceNumberForSlices; gdcm::Scanner::MappingType& tagValueMappings = const_cast(tagValueMappings_); //DICOM tags which should be added to the image properties const gdcm::Tag tagSliceLocation(0x0020, 0x1041); // slice location const gdcm::Tag tagInstanceNumber(0x0020, 0x0013); // (image) instance number const gdcm::Tag tagSOPInstanceNumber(0x0008, 0x0018); // SOP instance number unsigned int timeStep(0); std::string propertyKeySliceLocation = "dicom.image.0020.1041"; std::string propertyKeyInstanceNumber = "dicom.image.0020.0013"; std::string propertyKeySOPInstanceNumber = "dicom.image.0008.0018"; // tags for each image for ( std::list::iterator i = imageBlock.begin(); i != imageBlock.end(); i++, timeStep++ ) { const StringContainer& files = (*i); unsigned int slice(0); for ( StringContainer::const_iterator fIter = files.begin(); fIter != files.end(); ++fIter, ++slice ) { filesForSlices.SetTableValue( slice, *fIter ); gdcm::Scanner::TagToValue tagValueMapForFile = tagValueMappings[fIter->c_str()]; if(tagValueMapForFile.find(tagSliceLocation) != tagValueMapForFile.end()) sliceLocationForSlices.SetTableValue(slice, tagValueMapForFile[tagSliceLocation]); if(tagValueMapForFile.find(tagInstanceNumber) != tagValueMapForFile.end()) instanceNumberForSlices.SetTableValue(slice, tagValueMapForFile[tagInstanceNumber]); if(tagValueMapForFile.find(tagSOPInstanceNumber) != tagValueMapForFile.end()) SOPInstanceNumberForSlices.SetTableValue(slice, tagValueMapForFile[tagSOPInstanceNumber]); } image->SetProperty( "files", StringLookupTableProperty::New( filesForSlices ) ); //If more than one time step add postfix ".t" + timestep if(timeStep != 0) { std::ostringstream postfix; postfix << ".t" << timeStep; propertyKeySliceLocation.append(postfix.str()); propertyKeyInstanceNumber.append(postfix.str()); propertyKeySOPInstanceNumber.append(postfix.str()); } image->SetProperty( propertyKeySliceLocation.c_str(), StringLookupTableProperty::New( sliceLocationForSlices ) ); image->SetProperty( propertyKeyInstanceNumber.c_str(), StringLookupTableProperty::New( instanceNumberForSlices ) ); image->SetProperty( propertyKeySOPInstanceNumber.c_str(), StringLookupTableProperty::New( SOPInstanceNumberForSlices ) ); } // Copy tags for series, study, patient level (leave interpretation to application). // These properties will be copied to the DataNode by DicomSeriesReader. // tags for the series (we just use the one that ITK copied to its dictionary (proably that of the last slice) const itk::MetaDataDictionary& dict = io->GetMetaDataDictionary(); const TagToPropertyMapType& propertyLookup = DicomSeriesReader::GetDICOMTagsToMITKPropertyMap(); itk::MetaDataDictionary::ConstIterator dictIter = dict.Begin(); while ( dictIter != dict.End() ) { //MITK_DEBUG << "Key " << dictIter->first; std::string value; if ( itk::ExposeMetaData( dict, dictIter->first, value ) ) { //MITK_DEBUG << "Value " << value; TagToPropertyMapType::const_iterator valuePosition = propertyLookup.find( dictIter->first ); if ( valuePosition != propertyLookup.end() ) { std::string propertyKey = valuePosition->second; //MITK_DEBUG << "--> " << propertyKey; image->SetProperty( propertyKey.c_str(), StringProperty::New(value) ); } } else { MITK_WARN << "Tag " << dictIter->first << " not read as string as expected. Ignoring..." ; } ++dictIter; } // copy imageblockdescriptor as properties image->SetProperty("dicomseriesreader.SOPClass", StringProperty::New(blockInfo.GetSOPClassUIDAsString())); image->SetProperty("dicomseriesreader.ReaderImplementationLevelString", StringProperty::New(ReaderImplementationLevelToString( blockInfo.GetReaderImplementationLevel() ))); image->SetProperty("dicomseriesreader.ReaderImplementationLevel", GenericProperty::New( blockInfo.GetReaderImplementationLevel() )); image->SetProperty("dicomseriesreader.PixelSpacingInterpretationString", StringProperty::New(PixelSpacingInterpretationToString( blockInfo.GetPixelSpacingType() ))); image->SetProperty("dicomseriesreader.PixelSpacingInterpretation", GenericProperty::New(blockInfo.GetPixelSpacingType())); image->SetProperty("dicomseriesreader.MultiFrameImage", BoolProperty::New(blockInfo.IsMultiFrameImage())); image->SetProperty("dicomseriesreader.GantyTiltCorrected", BoolProperty::New(blockInfo.HasGantryTiltCorrected())); image->SetProperty("dicomseriesreader.3D+t", BoolProperty::New(blockInfo.HasMultipleTimePoints())); } void DicomSeriesReader::FixSpacingInformation( mitk::Image* image, const ImageBlockDescriptor& imageBlockDescriptor ) { // spacing provided by ITK/GDCM Vector3D imageSpacing = image->GetGeometry()->GetSpacing(); ScalarType imageSpacingX = imageSpacing[0]; ScalarType imageSpacingY = imageSpacing[1]; // spacing as desired by MITK (preference for "in patient", else "on detector", or "1.0/1.0") ScalarType desiredSpacingX = imageSpacingX; ScalarType desiredSpacingY = imageSpacingY; imageBlockDescriptor.GetDesiredMITKImagePixelSpacing( desiredSpacingX, desiredSpacingY ); MITK_DEBUG << "Loaded spacing: " << imageSpacingX << "/" << imageSpacingY; MITK_DEBUG << "Corrected spacing: " << desiredSpacingX << "/" << desiredSpacingY; imageSpacing[0] = desiredSpacingX; imageSpacing[1] = desiredSpacingY; image->GetGeometry()->SetSpacing( imageSpacing ); } -} // end namespace mitk +void DicomSeriesReader::LoadDicom(const StringContainer &filenames, DataNode &node, bool sort, bool load4D, bool correctTilt, UpdateCallBackMethod callback, Image::Pointer preLoadedImageBlock) +{ + const char* previousCLocale = setlocale(LC_NUMERIC, NULL); + setlocale(LC_NUMERIC, "C"); + std::locale previousCppLocale( std::cin.getloc() ); + std::locale l( "C" ); + std::cin.imbue(l); + + ImageBlockDescriptor imageBlockDescriptor; + + const gdcm::Tag tagImagePositionPatient(0x0020,0x0032); // Image Position (Patient) + const gdcm::Tag tagImageOrientation(0x0020, 0x0037); // Image Orientation + const gdcm::Tag tagSeriesInstanceUID(0x0020, 0x000e); // Series Instance UID + const gdcm::Tag tagSOPClassUID(0x0008, 0x0016); // SOP class UID + const gdcm::Tag tagModality(0x0008, 0x0060); // modality + const gdcm::Tag tagPixelSpacing(0x0028, 0x0030); // pixel spacing + const gdcm::Tag tagImagerPixelSpacing(0x0018, 0x1164); // imager pixel spacing + const gdcm::Tag tagNumberOfFrames(0x0028, 0x0008); // number of frames + + try + { + Image::Pointer image = preLoadedImageBlock.IsNull() ? Image::New() : preLoadedImageBlock; + CallbackCommand *command = callback ? new CallbackCommand(callback) : 0; + bool initialize_node = false; + + /* special case for Philips 3D+t ultrasound images */ + if ( DicomSeriesReader::IsPhilips3DDicom(filenames.front().c_str()) ) + { + // TODO what about imageBlockDescriptor? + // TODO what about preLoadedImageBlock? + ReadPhilips3DDicom(filenames.front().c_str(), image); + initialize_node = true; + } + else + { + /* default case: assume "normal" image blocks, possibly 3D+t */ + bool canLoadAs4D(true); + gdcm::Scanner scanner; + ScanForSliceInformation(filenames, scanner); + + // need non-const access for map + gdcm::Scanner::MappingType& tagValueMappings = const_cast(scanner.GetMappings()); + + std::list imageBlocks = SortIntoBlocksFor3DplusT( filenames, tagValueMappings, sort, canLoadAs4D ); + unsigned int volume_count = imageBlocks.size(); + + imageBlockDescriptor.SetSeriesInstanceUID( DicomSeriesReader::ConstCharStarToString( scanner.GetValue( filenames.front().c_str(), tagSeriesInstanceUID ) ) ); + imageBlockDescriptor.SetSOPClassUID( DicomSeriesReader::ConstCharStarToString( scanner.GetValue( filenames.front().c_str(), tagSOPClassUID ) ) ); + imageBlockDescriptor.SetModality( DicomSeriesReader::ConstCharStarToString( scanner.GetValue( filenames.front().c_str(), tagModality ) ) ); + imageBlockDescriptor.SetNumberOfFrames( ConstCharStarToString( scanner.GetValue( filenames.front().c_str(), tagNumberOfFrames ) ) ); + imageBlockDescriptor.SetPixelSpacingInformation( ConstCharStarToString( scanner.GetValue( filenames.front().c_str(), tagPixelSpacing ) ), + ConstCharStarToString( scanner.GetValue( filenames.front().c_str(), tagImagerPixelSpacing ) ) ); + + GantryTiltInformation tiltInfo; + + // check possibility of a single slice with many timesteps. In this case, don't check for tilt, no second slice possible + if ( !imageBlocks.empty() && imageBlocks.front().size() > 1 && correctTilt) + { + // check tiltedness here, potentially fixup ITK's loading result by shifting slice contents + // check first and last position slice from tags, make some calculations to detect tilt + + std::string firstFilename(imageBlocks.front().front()); + // calculate from first and last slice to minimize rounding errors + std::string secondFilename(imageBlocks.front().back()); + + std::string imagePosition1( ConstCharStarToString( tagValueMappings[ firstFilename.c_str() ][ tagImagePositionPatient ] ) ); + std::string imageOrientation( ConstCharStarToString( tagValueMappings[ firstFilename.c_str() ][ tagImageOrientation ] ) ); + std::string imagePosition2( ConstCharStarToString( tagValueMappings[secondFilename.c_str() ][ tagImagePositionPatient ] ) ); + + bool ignoredConversionError(-42); // hard to get here, no graceful way to react + Point3D origin1( DICOMStringToPoint3D( imagePosition1, ignoredConversionError ) ); + Point3D origin2( DICOMStringToPoint3D( imagePosition2, ignoredConversionError ) ); + + Vector3D right; right.Fill(0.0); + Vector3D up; right.Fill(0.0); // might be down as well, but it is just a name at this point + DICOMStringToOrientationVectors( imageOrientation, right, up, ignoredConversionError ); + + tiltInfo = GantryTiltInformation ( origin1, origin2, right, up, filenames.size()-1 ); + correctTilt = tiltInfo.IsSheared() && tiltInfo.IsRegularGantryTilt(); + } + else + { + correctTilt = false; // we CANNOT do that + } + + imageBlockDescriptor.SetHasGantryTiltCorrected( correctTilt ); + + if (volume_count == 1 || !canLoadAs4D || !load4D) + { + + DcmIoType::Pointer io; + image = MultiplexLoadDICOMByITK( imageBlocks.front(), correctTilt, tiltInfo, io, command, preLoadedImageBlock ); // load first 3D block + + imageBlockDescriptor.AddFiles(imageBlocks.front()); // only the first part is loaded + imageBlockDescriptor.SetHasMultipleTimePoints( false ); -#include + FixSpacingInformation( image, imageBlockDescriptor ); + CopyMetaDataToImageProperties( imageBlocks.front(), scanner.GetMappings(), io, imageBlockDescriptor, image); + + initialize_node = true; + } + else if (volume_count > 1) + { + imageBlockDescriptor.AddFiles(filenames); // all is loaded + imageBlockDescriptor.SetHasMultipleTimePoints( true ); + + DcmIoType::Pointer io; + image = MultiplexLoadDICOMByITK4D( imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command, preLoadedImageBlock ); + + initialize_node = true; + } + } + + if (initialize_node) + { + // forward some image properties to node + node.GetPropertyList()->ConcatenatePropertyList( image->GetPropertyList(), true ); + + node.SetData( image ); + setlocale(LC_NUMERIC, previousCLocale); + std::cin.imbue(previousCppLocale); + } + + MITK_DEBUG << "--------------------------------------------------------------------------------"; + MITK_DEBUG << "DICOM files loaded (from series UID " << imageBlockDescriptor.GetSeriesInstanceUID() << "):"; + MITK_DEBUG << " " << imageBlockDescriptor.GetFilenames().size() << " '" << imageBlockDescriptor.GetModality() << "' files (" << imageBlockDescriptor.GetSOPClassUIDAsString() << ") loaded into 1 mitk::Image"; + MITK_DEBUG << " multi-frame: " << (imageBlockDescriptor.IsMultiFrameImage()?"Yes":"No"); + MITK_DEBUG << " reader support: " << ReaderImplementationLevelToString(imageBlockDescriptor.GetReaderImplementationLevel()); + MITK_DEBUG << " pixel spacing type: " << PixelSpacingInterpretationToString( imageBlockDescriptor.GetPixelSpacingType() ) << " " << image->GetGeometry()->GetSpacing()[0] << "/" << image->GetGeometry()->GetSpacing()[0]; + MITK_DEBUG << " gantry tilt corrected: " << (imageBlockDescriptor.HasGantryTiltCorrected()?"Yes":"No"); + MITK_DEBUG << " 3D+t: " << (imageBlockDescriptor.HasMultipleTimePoints()?"Yes":"No"); + MITK_DEBUG << "--------------------------------------------------------------------------------"; + } + catch (std::exception& e) + { + // reset locale then throw up + setlocale(LC_NUMERIC, previousCLocale); + std::cin.imbue(previousCppLocale); + + MITK_DEBUG << "Caught exception in DicomSeriesReader::LoadDicom"; + + throw e; + } +} + +void +DicomSeriesReader::ScanForSliceInformation(const StringContainer &filenames, gdcm::Scanner& scanner) +{ + const gdcm::Tag tagImagePositionPatient(0x0020,0x0032); //Image position (Patient) + scanner.AddTag(tagImagePositionPatient); + + const gdcm::Tag tagSeriesInstanceUID(0x0020, 0x000e); // Series Instance UID + scanner.AddTag(tagSeriesInstanceUID); + + const gdcm::Tag tagImageOrientation(0x0020,0x0037); //Image orientation + scanner.AddTag(tagImageOrientation); + + const gdcm::Tag tagSliceLocation(0x0020, 0x1041); // slice location + scanner.AddTag( tagSliceLocation ); + + const gdcm::Tag tagInstanceNumber(0x0020, 0x0013); // (image) instance number + scanner.AddTag( tagInstanceNumber ); + + const gdcm::Tag tagSOPInstanceNumber(0x0008, 0x0018); // SOP instance number + scanner.AddTag( tagSOPInstanceNumber ); + + const gdcm::Tag tagPixelSpacing(0x0028, 0x0030); // Pixel Spacing + scanner.AddTag( tagPixelSpacing ); + + const gdcm::Tag tagImagerPixelSpacing(0x0018, 0x1164); // Imager Pixel Spacing + scanner.AddTag( tagImagerPixelSpacing ); + + const gdcm::Tag tagModality(0x0008, 0x0060); // Modality + scanner.AddTag( tagModality ); + + const gdcm::Tag tagSOPClassUID(0x0008, 0x0016); // SOP Class UID + scanner.AddTag( tagSOPClassUID ); + + const gdcm::Tag tagNumberOfFrames(0x0028, 0x0008); // number of frames + scanner.AddTag( tagNumberOfFrames ); + + scanner.Scan(filenames); // make available image information for each file +} + +std::list +DicomSeriesReader::SortIntoBlocksFor3DplusT( + const StringContainer& presortedFilenames, + const gdcm::Scanner::MappingType& tagValueMappings, + bool /*sort*/, + bool& canLoadAs4D ) +{ + std::list imageBlocks; + + // ignore sort request, because most likely re-sorting is now needed due to changes in GetSeries(bug #8022) + StringContainer sorted_filenames = DicomSeriesReader::SortSeriesSlices(presortedFilenames); + + std::string firstPosition; + unsigned int numberOfBlocks(0); // number of 3D image blocks + + static const gdcm::Tag tagImagePositionPatient(0x0020,0x0032); //Image position (Patient) + + // loop files to determine number of image blocks + for (StringContainer::const_iterator fileIter = sorted_filenames.begin(); + fileIter != sorted_filenames.end(); + ++fileIter) + { + gdcm::Scanner::TagToValue tagToValueMap = tagValueMappings.find( fileIter->c_str() )->second; + + if(tagToValueMap.find(tagImagePositionPatient) == tagToValueMap.end()) + { + // we expect to get images w/ missing position information ONLY as separated blocks. + assert( presortedFilenames.size() == 1 ); + numberOfBlocks = 1; + break; + } + + std::string position = tagToValueMap.find(tagImagePositionPatient)->second; + MITK_DEBUG << " " << *fileIter << " at " << position; + if (firstPosition.empty()) + { + firstPosition = position; + } + + if ( position == firstPosition ) + { + ++numberOfBlocks; + } + else + { + break; // enough information to know the number of image blocks + } + } + + MITK_DEBUG << " ==> Assuming " << numberOfBlocks << " time steps"; + + if (numberOfBlocks == 0) return imageBlocks; // only possible if called with no files + + + // loop files to sort them into image blocks + unsigned int numberOfExpectedSlices(0); + for (unsigned int block = 0; block < numberOfBlocks; ++block) + { + StringContainer filesOfCurrentBlock; + + for ( StringContainer::const_iterator fileIter = sorted_filenames.begin() + block; + fileIter != sorted_filenames.end(); + //fileIter += numberOfBlocks) // TODO shouldn't this work? give invalid iterators on first attempts + ) + { + filesOfCurrentBlock.push_back( *fileIter ); + for (unsigned int b = 0; b < numberOfBlocks; ++b) + { + if (fileIter != sorted_filenames.end()) + ++fileIter; + } + } + + imageBlocks.push_back(filesOfCurrentBlock); + + if (block == 0) + { + numberOfExpectedSlices = filesOfCurrentBlock.size(); + } + else + { + if (filesOfCurrentBlock.size() != numberOfExpectedSlices) + { + MITK_WARN << "DicomSeriesReader expected " << numberOfBlocks + << " image blocks of " + << numberOfExpectedSlices + << " images each. Block " + << block + << " got " + << filesOfCurrentBlock.size() + << " instead. Cannot load this as 3D+t"; // TODO implement recovery (load as many slices 3D+t as much as possible) + canLoadAs4D = false; + } + } + } + + return imageBlocks; +} + +Image::Pointer +DicomSeriesReader +::MultiplexLoadDICOMByITK(const StringContainer& filenames, bool correctTilt, const GantryTiltInformation& tiltInfo, DcmIoType::Pointer& io, CallbackCommand* command, Image::Pointer preLoadedImageBlock) +{ + io = DcmIoType::New(); + io->SetFileName(filenames.front().c_str()); + io->ReadImageInformation(); + + if (io->GetPixelType() == itk::ImageIOBase::SCALAR) + { + return MultiplexLoadDICOMByITKScalar(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + } + else if (io->GetPixelType() == itk::ImageIOBase::RGB) + { + return MultiplexLoadDICOMByITKRGBPixel(filenames, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + } + else + { + return NULL; + } +} + +Image::Pointer +DicomSeriesReader +::MultiplexLoadDICOMByITK4D( std::list& imageBlocks, ImageBlockDescriptor imageBlockDescriptor, bool correctTilt, const GantryTiltInformation& tiltInfo, DcmIoType::Pointer& io, CallbackCommand* command, Image::Pointer preLoadedImageBlock) +{ + io = DcmIoType::New(); + io->SetFileName(imageBlocks.front().front().c_str()); + io->ReadImageInformation(); + + if (io->GetPixelType() == itk::ImageIOBase::SCALAR) + { + return MultiplexLoadDICOMByITK4DScalar(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + } + else if (io->GetPixelType() == itk::ImageIOBase::RGB) + { + return MultiplexLoadDICOMByITK4DRGBPixel(imageBlocks, imageBlockDescriptor, correctTilt, tiltInfo, io, command ,preLoadedImageBlock); + } + else + { + return NULL; + } +} + +} // end namespace mitk diff --git a/Core/Code/IO/mitkDicomSeriesReader.h b/Core/Code/IO/mitkDicomSeriesReader.h index 7857e3e04d..ca6bbc2c2d 100644 --- a/Core/Code/IO/mitkDicomSeriesReader.h +++ b/Core/Code/IO/mitkDicomSeriesReader.h @@ -1,939 +1,967 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef mitkDicomSeriesReader_h #define mitkDicomSeriesReader_h #include "mitkDataNode.h" #include "mitkConfig.h" #include #include #include #ifdef NOMINMAX # define DEF_NOMINMAX # undef NOMINMAX #endif #include #ifdef DEF_NOMINMAX # ifndef NOMINMAX # define NOMINMAX # endif # undef DEF_NOMINMAX #endif #include #include namespace mitk { /** \brief Loading DICOM images as MITK images. - \ref DicomSeriesReader_purpose - \ref DicomSeriesReader_limitations - \ref DicomSeriesReader_usage - \ref DicomSeriesReader_sorting - \ref DicomSeriesReader_sorting1 - \ref DicomSeriesReader_sorting2 - \ref DicomSeriesReader_sorting3 - \ref DicomSeriesReader_sorting4 - \ref DicomSeriesReader_gantrytilt - \ref DicomSeriesReader_pixelspacing - \ref DicomSeriesReader_nextworkitems - \ref DicomSeriesReader_whynotinitk - \ref DicomSeriesReader_tests \section DicomSeriesReader_purpose Purpose DicomSeriesReader serves as a central class for loading DICOM images as mitk::Image. As the term "DICOM image" covers a huge variety of possible modalities and implementations, and since MITK assumes that 3D images are made up of continuous blocks of slices without any gaps or changes in orientation, the loading mechanism must implement a number of decisions and compromises. The main intention of this implementation is not efficiency but correctness of generated slice positions and pixel spacings! \section DicomSeriesReader_limitations Assumptions and limitations The class is working only with GDCM 2.0.14 (or possibly newer). This version is the default of an MITK super-build. Support for other versions or ITK's DicomIO was dropped because of the associated complexity of DicomSeriesReader. \b Assumptions - expected to work with certain SOP Classes (mostly CT Image Storage and MR Image Storage) - see ImageBlockDescriptor.GetReaderImplementationLevel() method for the details - special treatment for a certain type of Philips 3D ultrasound (recogized by tag 3001,0010 set to "Philips3D") - loader will always attempt to read multiple single slices as a single 3D image volume (i.e. mitk::Image) - slices will be grouped by basic properties such as orientation, rows, columns, spacing and grouped into as large blocks as possible - images which do NOT report a position or orientation in space (Image Position Patient, Image Orientation) will be assigned defaults - image position (0,0,0) - image orientation (1,0,0), (0,1,0) - such images will always be grouped separately since spatial grouping / sorting makes no sense for them \b Options - images that cover the same piece of space (i.e. position, orientation, and dimensions are equal) can be interpreted as time-steps of the same image, i.e. a series will be loaded as 3D+t \b Limitations - the 3D+t assumption only works if all time-steps have an equal number of slices and if all have the Acquisition Time attribute set to meaningful values \section DicomSeriesReader_usage Usage The starting point for an application is a set of DICOM files that should be loaded. For convenience, DicomSeriesReader can also parse a whole directory for DICOM files, but an application should better know exactly what to load. Loading is then done in two steps: 1. Group the files into spatial blocks by calling GetSeries(). This method will sort all passed files into meaningful blocks that could fit into an mitk::Image. Sorting for 3D+t loading is optional but default. The \b return value of this function is a list of descriptors, which describe a grouped list of files with its most basic properties: - SOP Class (CT Image Storage, Secondary Capture Image Storage, etc.) - Modality - What type of pixel spacing can be read from the provided DICOM tags - How well DicomSeriesReader is prepared to load this type of data 2. Load a sorted set of files by calling LoadDicomSeries(). This method expects go receive the sorting output of GetSeries(). The method will then invoke ITK methods configured with GDCM-IO classes to actually load the files into memory and put them into mitk::Images. Again, loading as 3D+t is optional. Example: \code // only a directory is known at this point: /home/who/dicom DicomSeriesReader::FileNamesGrouping allImageBlocks = DicomSeriesReader::GetSeries("/home/who/dicom/"); // file now divided into groups of identical image size, orientation, spacing, etc. // each of these lists should be loadable as an mitk::Image. DicomSeriesReader::StringContainer seriesToLoad = allImageBlocks[...]; // decide what to load // final step: load into DataNode (can result in 3D+t image) DataNode::Pointer node = DicomSeriesReader::LoadDicomSeries( oneBlockSorted ); Image::Pointer image = dynamic_cast( node->GetData() ); \endcode \section DicomSeriesReader_sorting Logic for sorting 2D slices from DICOM images into 3D+t blocks for mitk::Image The general sorting mechanism (implemented in GetSeries) groups and sorts a set of DICOM files, each assumed to contain a single CT/MR slice. In the following we refer to those file groups as "blocks", since this is what they are meant to become when loaded into an mitk::Image. \subsection DicomSeriesReader_sorting1 Step 1: Avoiding pure non-sense A first pass separates slices that cannot possibly be loaded together because of restrictions of mitk::Image. After this steps, each block contains only slices that match in all of the following DICOM tags: - (0020,000e) Series Instance UID - (0020,0037) Image Orientation - (0028,0030) Pixel Spacing - (0018,1164) Imager Pixel Spacing - (0018,0050) Slice Thickness - (0028,0010) Number Of Rows - (0028,0011) Number Of Columns - (0028,0008) Number Of Frames \subsection DicomSeriesReader_sorting2 Step 2: Sort slices spatially Before slices are further analyzed, they are sorted spatially. As implemented by GdcmSortFunction(), slices are sorted by 1. distance from origin (calculated using (0020,0032) Image Position Patient and (0020,0037) Image Orientation) 2. when distance is equal, (0020,0012) Aquisition Number, (0008,0032) Acquisition Time and (0018,1060) Trigger Time are used as a backup criterions (necessary for meaningful 3D+t sorting) \subsection DicomSeriesReader_sorting3 Step 3: Ensure equal z spacing Since inter-slice distance is not recorded in DICOM tags, we must ensure that blocks are made up of slices that have equal distances between neighboring slices. This is especially necessary because itk::ImageSeriesReader is later used for the actual loading, and this class expects (and does nocht verify) equal inter-slice distance (see \ref DicomSeriesReader_whatweknowaboutitk). To achieve such grouping, the inter-slice distance is calculated from the first two different slice positions of a block. Following slices are added to a block as long as they can be added by adding the calculated inter-slice distance to the last slice of the block. Slices that do not fit into the expected distance pattern, are set aside for further analysis. This grouping is done until each file has been assigned to a group. Slices that share a position in space are also sorted into separate blocks during this step. So the result of this step is a set of blocks that contain only slices with equal z spacing and uniqe slices at each position. \subsection DicomSeriesReader_sorting4 Step 4 (optional): group 3D blocks as 3D+t when possible This last step depends on an option of GetSeries(). When requested, image blocks from the previous step are merged again whenever two blocks occupy the same portion of space (i.e. same origin, number of slices and z-spacing). \section DicomSeriesReader_gantrytilt Handling of gantry tilt When CT gantry tilt is used, the gantry plane (= X-Ray source and detector ring) and the vertical plane do not align anymore. This scanner feature is used for example to reduce metal artifacs (e.g. Lee C , Evaluation of Using CT Gantry Tilt Scan on Head and Neck Cancer Patients with Dental Structure: Scans Show Less Metal Artifacts. Presented at: Radiological Society of North America 2011 Scientific Assembly and Annual Meeting; November 27- December 2, 2011 Chicago IL.). The acquired planes of such CT series do not match the expectations of a orthogonal geometry in mitk::Image: if you stack the slices, they show a small shift along the Y axis: \verbatim without tilt with tilt |||||| ////// |||||| ////// -- |||||| --------- ////// -------- table orientation |||||| ////// |||||| ////// Stacked slices: without tilt with tilt -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- -------------- \endverbatim As such gemetries do not in conjunction with mitk::Image, DicomSeriesReader performs a correction for such series if the groupImagesWithGantryTilt or correctGantryTilt flag in GetSeries and LoadDicomSeries is set (default = on). The correction algorithms undoes two errors introduced by ITK's ImageSeriesReader: - the plane shift that is ignored by ITK's reader is recreated by applying a shearing transformation using itk::ResampleFilter. - the spacing is corrected (it is calculated by ITK's reader from the distance between two origins, which is NOT the slice distance in this special case) Both errors are introduced in itkImageSeriesReader.txx (ImageSeriesReader::GenerateOutputInformation(void)), lines 176 to 245 (as of ITK 3.20) For the correction, we examine two consecutive slices of a series, both described as a pair (origin/orientation): - we calculate if the first origin is on a line along the normal of the second slice - if this is not the case, the geometry will not fit a normal mitk::Image/mitk::Geometry3D - we then project the second origin into the first slice's coordinate system to quantify the shift - both is done in class GantryTiltInformation with quite some comments. The geometry of image stacks with tilted geometries is illustrated below: - green: the DICOM images as described by their tags: origin as a point with the line indicating the orientation - red: the output of ITK ImageSeriesReader: wrong, larger spacing, no tilt - blue: how much a shear must correct \image tilt-correction.jpg \section DicomSeriesReader_whatweknowaboutitk The actual image loading process When calling LoadDicomSeries(), this method "mainly" uses an instance of itk::ImageSeriesReader, configured with an itk::GDCMImageIO object. Because DicomSeriesReader works around some of the behaviors of these classes, the following is a list of features that we find in the code and need to work with: - itk::ImageSeriesReader::GenerateOutputInformation() does the z-spacing handling + spacing is directly determined by comparing (euclidean distance) the origins of the first two slices of a series * this is GOOD because there is no reliable z-spacing information in DICOM images * this is bad because it does not work with gantry tilt, in which case the slice distance is SMALLER than the distance between two origins (see section on tilt) - origin and spacing are calculated by GDCMImageIO and re-used in itk::ImageSeriesReader + the origins are read from appropriate tags, nothing special about that + the spacing is read by gdcm::ImageReader, gdcm::ImageHelper::GetSpacingValue() from a tag determined by gdcm::ImageHelper::GetSpacingTagFromMediaStorage(), which basically determines ONE appropriate pixel spacing tag for each media storage type (ct image, mr image, secondary capture image, etc.) * this is fine for modalities such as CT/MR where the "Pixel Spacing" tag is mandatory, but for other modalities such as CR or Secondary Capture, the tag "Imager Pixel Spacing" is taken, which is no only optional but also has a more complicated relation with the "Pixel Spacing" tag. For this reason we check/modify the pixel spacing reported by itk::ImageSeriesReader after loading the image (see \ref DicomSeriesReader_pixelspacing) AFTER loading, DicomSeriesReader marks some of its findings as mitk::Properties to the loaded Image and DataNode: - dicomseriesreader.SOPClass : DICOM SOP Class as readable string (instead of a UID) - dicomseriesreader.ReaderImplementationLevelString : Confidence /Support level of the reader for this image as readable string - dicomseriesreader.ReaderImplementationLevel : Confidence /Support level of the reader for this image as enum value of type ReaderImplementationLevel - dicomseriesreader.PixelSpacingInterpretationString : Appropriate interpreteation of pixel spacing for this Image as readable string - dicomseriesreader.PixelSpacingInterpretation : Appropriate interpreteation of pixel spacing for this Image as enum value of type PixelSpacingInterpretation - dicomseriesreader.MultiFrameImage : bool flag to mark multi-frame images - dicomseriesreader.GantyTiltCorrected : bool flag to mark images where a gantry tilt was corrected to fit slices into an mitk::Image - dicomseriesreader.3D+t : bool flag to mark images with a time dimension (multiple 3D blocks of the same size at the same position in space) \section DicomSeriesReader_pixelspacing Handling of pixel spacing The reader implementes what is described in DICOM Part 3, chapter 10.7 (Basic Pixel Spacing Calibration Macro): Both tags - (0028,0030) Pixel Spacing and - (0018,1164) Imager Pixel Spacing are evaluated and the pixel spacing is set to the spacing within the patient when tags allow that. The result of pixel spacing interpretation can be read from a property "dicomseriesreader.PixelSpacingInterpretation", which refers to one of the enumerated values of type PixelSpacingInterpretation; \section DicomSeriesReader_supportedmodalities Limitations for specific modalities - Enhanced Computed Tomography / Magnetic Resonance Images are currently NOT supported at all, because we lack general support for multi-frame images. - Nuclear Medicine Images are not supported fully supported, only the single-frame variants are loaded properly. \section DicomSeriesReader_nextworkitems Possible enhancements This is a short list of ideas for enhancement: - Class has historically grown and should be reviewed again. There is probably too many duplicated scanning code - Multi-frame images don't mix well with the curent assumption of "one file - one slice", which is assumed by our code - It should be checked how well GDCM and ITK support these files (some load, some don't) - Specializations such as the Philips 3D code should be handled in a more generic way. The current handling of Philips 3D images is not nice at all \section DicomSeriesReader_whynotinitk Why is this not in ITK? Some of this code would probably be better located in ITK. It is just a matter of resources that this is not the case yet. Any attempts into this direction are welcome and can be supported. At least the gantry tilt correction should be a simple addition to itk::ImageSeriesReader. \section DicomSeriesReader_tests Tests regarding DICOM loading A number of tests have been implemented to check our assumptions regarding DICOM loading. Please see \ref DICOMTesting \todo refactor all the protected helper objects/methods into a separate header so we compile faster */ class MITK_CORE_EXPORT DicomSeriesReader { public: /** \brief Lists of filenames. */ typedef std::vector StringContainer; /** \brief Interface for the progress callback. */ typedef void (*UpdateCallBackMethod)(float); /** \brief Describes how well the reader is tested for a certain file type. Applications should not rely on the outcome for images which are reported ReaderImplementationLevel_Implemented or ReaderImplementationLevel_Unsupported. Errors to load images which are reported as ReaderImplementationLevel_Supported are considered bugs. For ReaderImplementationLevel_PartlySupported please check the appropriate paragraph in \ref DicomSeriesReader_supportedmodalities */ typedef enum { ReaderImplementationLevel_Supported, /// loader code and tests are established ReaderImplementationLevel_PartlySupported, /// loader code and tests are establised for specific parts of a SOP Class ReaderImplementationLevel_Implemented, /// loader code is implemented but not accompanied by tests ReaderImplementationLevel_Unsupported, /// loader code is not working with this SOP Class } ReaderImplementationLevel; /** \brief How the mitk::Image spacing should be interpreted. Compare DICOM PS 3.3 10.7 (Basic Pixel Spacing Calibration Macro). */ typedef enum { PixelSpacingInterpretation_SpacingInPatient, /// distances are mm within a patient PixelSpacingInterpretation_SpacingAtDetector, /// distances are mm at detector surface PixelSpacingInterpretation_SpacingUnknown /// NO spacing information is present, we use (1,1) as default } PixelSpacingInterpretation; /** \brief Return type of GetSeries, describes a logical group of files. Files grouped into a single 3D or 3D+t block are described by an instance of this class. Relevant descriptive properties can be used to provide the application user with meaningful choices. */ class MITK_CORE_EXPORT ImageBlockDescriptor { public: /// List of files in this group StringContainer GetFilenames() const; /// A unique ID describing this bloc (enhanced Series Instance UID). std::string GetImageBlockUID() const; /// The Series Instance UID. std::string GetSeriesInstanceUID() const; /// Series Modality (CT, MR, etc.) std::string GetModality() const; /// SOP Class UID as readable string (Computed Tomography Image Storage, Secondary Capture Image Storage, etc.) std::string GetSOPClassUIDAsString() const; /// SOP Class UID as DICOM UID std::string GetSOPClassUID() const; /// Confidence of the reader that this block can be read successfully. ReaderImplementationLevel GetReaderImplementationLevel() const; /// Whether or not the block contains a gantry tilt which will be "corrected" during loading bool HasGantryTiltCorrected() const; /// Whether or not mitk::Image spacing relates to the patient bool PixelSpacingRelatesToPatient() const; /// Whether or not mitk::Image spacing relates to the detector surface bool PixelSpacingRelatesToDetector() const; /// Whether or not mitk::Image spacing is of unknown origin bool PixelSpacingIsUnknown() const; /// How the mitk::Image spacing can meaningfully be interpreted. PixelSpacingInterpretation GetPixelSpacingType() const; /// 3D+t or not bool HasMultipleTimePoints() const; /// Multi-frame image(s) or not bool IsMultiFrameImage() const; ImageBlockDescriptor(); ~ImageBlockDescriptor(); private: friend class DicomSeriesReader; ImageBlockDescriptor(const StringContainer& files); void AddFile(const std::string& file); void AddFiles(const StringContainer& files); void SetImageBlockUID(const std::string& uid); void SetSeriesInstanceUID(const std::string& uid); void SetModality(const std::string& modality); void SetNumberOfFrames(const std::string& ); void SetSOPClassUID(const std::string& mediaStorageSOPClassUID); void SetHasGantryTiltCorrected(bool); void SetPixelSpacingInformation(const std::string& pixelSpacing, const std::string& imagerPixelSpacing); void SetHasMultipleTimePoints(bool); void GetDesiredMITKImagePixelSpacing( float& spacingX, float& spacingY) const; StringContainer m_Filenames; std::string m_ImageBlockUID; std::string m_SeriesInstanceUID; std::string m_Modality; std::string m_SOPClassUID; bool m_HasGantryTiltCorrected; std::string m_PixelSpacing; std::string m_ImagerPixelSpacing; bool m_HasMultipleTimePoints; bool m_IsMultiFrameImage; }; typedef std::map FileNamesGrouping; /** \brief Provide combination of preprocessor defines that was active during compilation. Since this class is a combination of several possible implementations, separated only by ifdef's, calling instances might want to know which flags were active at compile time. */ static std::string GetConfigurationString(); /** \brief Checks if a specific file contains DICOM data. */ static bool IsDicom(const std::string &filename); /** \brief see other GetSeries(). Find all series (and sub-series -- see details) in a particular directory. */ static FileNamesGrouping GetSeries(const std::string &dir, bool groupImagesWithGantryTilt, const StringContainer &restrictions = StringContainer()); /** \brief see other GetSeries(). \warning Untested, could or could not work. This differs only by having an additional restriction to a single known DICOM series. Internally, it uses the other GetSeries() method. */ static StringContainer GetSeries(const std::string &dir, const std::string &series_uid, bool groupImagesWithGantryTilt, const StringContainer &restrictions = StringContainer()); /** \brief PREFERRED version of this method - scan and sort DICOM files. Parse a list of files for images of DICOM series. For each series, an enumeration of the files contained in it is created. \return The resulting maps UID-like keys (based on Series Instance UID and slice properties) to sorted lists of file names. SeriesInstanceUID will be enhanced to be unique for each set of file names that is later loadable as a single mitk::Image. This implies that Image orientation, slice thickness, pixel spacing, rows, and columns must be the same for each file (i.e. the image slice contained in the file). If this separation logic requires that a SeriesInstanceUID must be made more specialized, it will follow the same logic as itk::GDCMSeriesFileNames to enhance the UID with more digits and dots. Optionally, more tags can be used to separate files into different logical series by setting the restrictions parameter. \warning Adding restrictions is not yet implemented! */ static FileNamesGrouping GetSeries(const StringContainer& files, bool sortTo3DPlust, bool groupImagesWithGantryTilt, const StringContainer &restrictions = StringContainer()); /** \brief See other GetSeries(). Use GetSeries(const StringContainer& files, bool sortTo3DPlust, const StringContainer &restrictions) instead. */ static FileNamesGrouping GetSeries(const StringContainer& files, bool groupImagesWithGantryTilt, const StringContainer &restrictions = StringContainer()); /** Loads a DICOM series composed by the file names enumerated in the file names container. If a callback method is supplied, it will be called after every progress update with a progress value in [0,1]. \param filenames The filenames to load. \param sort Whether files should be sorted spatially (true) or not (false - maybe useful if presorted) \param load4D Whether to load the files as 3D+t (if possible) */ static DataNode::Pointer LoadDicomSeries(const StringContainer &filenames, bool sort = true, bool load4D = true, bool correctGantryTilt = true, UpdateCallBackMethod callback = 0, Image::Pointer preLoadedImageBlock = 0); /** \brief See LoadDicomSeries! Just a slightly different interface. If \p preLoadedImageBlock is provided, the reader will only "fake" loading and create appropriate mitk::Properties. */ static bool LoadDicomSeries(const StringContainer &filenames, DataNode &node, bool sort = true, bool load4D = true, bool correctGantryTilt = true, UpdateCallBackMethod callback = 0, Image::Pointer preLoadedImageBlock = 0); protected: /** \brief Return type of DicomSeriesReader::AnalyzeFileForITKImageSeriesReaderSpacingAssumption. Class contains the grouping result of method DicomSeriesReader::AnalyzeFileForITKImageSeriesReaderSpacingAssumption, which takes as input a number of images, which are all equally oriented and spatially sorted along their normal direction. The result contains of two blocks: a first one is the grouping result, all of those images can be loaded into one image block because they have an equal origin-to-origin distance without any gaps in-between. */ class SliceGroupingAnalysisResult { public: SliceGroupingAnalysisResult(); /** \brief Grouping result, all same origin-to-origin distance w/o gaps. */ StringContainer GetBlockFilenames(); /** \brief Remaining files, which could not be grouped. */ StringContainer GetUnsortedFilenames(); /** \brief Wheter or not the grouped result contain a gantry tilt. */ bool ContainsGantryTilt(); /** \brief Meant for internal use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption only. */ void AddFileToSortedBlock(const std::string& filename); /** \brief Meant for internal use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption only. */ void AddFileToUnsortedBlock(const std::string& filename); void AddFilesToUnsortedBlock(const StringContainer& filenames); /** \brief Meant for internal use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption only. \todo Could make sense to enhance this with an instance of GantryTiltInformation to store the whole result! */ void FlagGantryTilt(); /** \brief Only meaningful for use by AnalyzeFileForITKImageSeriesReaderSpacingAssumption. */ void UndoPrematureGrouping(); protected: StringContainer m_GroupedFiles; StringContainer m_UnsortedFiles; bool m_GantryTilt; }; /** \brief Gantry tilt analysis result. Takes geometry information for two slices of a DICOM series and calculates whether these fit into an orthogonal block or not. If NOT, they can either be the result of an acquisition with gantry tilt OR completly broken by some shearing transformation. Most calculations are done in the constructor, results can then be read via the remaining methods. */ class GantryTiltInformation { public: // two types to avoid any rounding errors typedef itk::Point Point3Dd; typedef itk::Vector Vector3Dd; /** \brief Just so we can create empty instances for assigning results later. */ GantryTiltInformation(); /** \brief THE constructor, which does all the calculations. Determining the amount of tilt is done by checking the distances of origin1 from planes through origin2. Two planes are considered: - normal vector along normal of slices (right x up): gives the slice distance - normal vector along orientation vector "up": gives the shift parallel to the plane orientation The tilt angle can then be calculated from these distances \param origin1 origin of the first slice \param origin2 origin of the second slice \param right right/up describe the orientatation of borth slices \param up right/up describe the orientatation of borth slices \param numberOfSlicesApart how many slices are the given origins apart (1 for neighboring slices) */ GantryTiltInformation( const Point3D& origin1, const Point3D& origin2, const Vector3D& right, const Vector3D& up, unsigned int numberOfSlicesApart); /** \brief Whether the slices were sheared. True if any of the shifts along right or up vector are non-zero. */ bool IsSheared() const; /** \brief Whether the shearing is a gantry tilt or more complicated. Gantry tilt will only produce shifts in ONE orientation, not in both. Since the correction code currently only coveres one tilt direction AND we don't know of medical images with two tilt directions, the loading code wants to check if our assumptions are true. */ bool IsRegularGantryTilt() const; /** \brief The offset distance in Y direction for each slice in mm (describes the tilt result). */ double GetMatrixCoefficientForCorrectionInWorldCoordinates() const; /** \brief The z / inter-slice spacing. Needed to correct ImageSeriesReader's result. */ double GetRealZSpacing() const; /** \brief The shift between first and last slice in mm. Needed to resize an orthogonal image volume. */ double GetTiltCorrectedAdditionalSize() const; /** \brief Calculated tilt angle in degrees. */ double GetTiltAngleInDegrees() const; protected: /** \brief Projection of point p onto line through lineOrigin in direction of lineDirection. */ Point3D projectPointOnLine( Point3Dd p, Point3Dd lineOrigin, Vector3Dd lineDirection ); double m_ShiftUp; double m_ShiftRight; double m_ShiftNormal; double m_ITKAssumedSliceSpacing; unsigned int m_NumberOfSlicesApart; }; /** \brief for internal sorting. */ typedef std::pair TwoStringContainers; /** \brief Maps DICOM tags to MITK properties. */ typedef std::map TagToPropertyMapType; /** \brief Ensure an equal z-spacing for a group of files. Takes as input a number of images, which are all equally oriented and spatially sorted along their normal direction. Internally used by GetSeries. Returns two lists: the first one contins slices of equal inter-slice spacing. The second list contains remaining files, which need to be run through AnalyzeFileForITKImageSeriesReaderSpacingAssumption again. Relevant code that is matched here is in itkImageSeriesReader.txx (ImageSeriesReader::GenerateOutputInformation(void)), lines 176 to 245 (as of ITK 3.20) */ static SliceGroupingAnalysisResult AnalyzeFileForITKImageSeriesReaderSpacingAssumption(const StringContainer& files, bool groupsOfSimilarImages, const gdcm::Scanner::MappingType& tagValueMappings_); /** \brief Safely convert const char* to std::string. */ static std::string ConstCharStarToString(const char* s); /** \brief Safely convert a string into pixel spacing x and y. */ static bool DICOMStringToSpacing(const std::string& s, float& spacingX, float& spacingY); /** \brief Convert DICOM string describing a point to Point3D. DICOM tags like ImagePositionPatient contain a position as float numbers separated by backslashes: \verbatim 42.7131\13.77\0.7 \endverbatim */ static Point3D DICOMStringToPoint3D(const std::string& s, bool& successful); /** \brief Convert DICOM string describing a point two Vector3D. DICOM tags like ImageOrientationPatient contain two vectors as float numbers separated by backslashes: \verbatim 42.7131\13.77\0.7\137.76\0.3 \endverbatim */ static void DICOMStringToOrientationVectors(const std::string& s, Vector3D& right, Vector3D& up, bool& successful); template static typename ImageType::Pointer // TODO this is NOT inplace! InPlaceFixUpTiltedGeometry( ImageType* input, const GantryTiltInformation& tiltInfo ); /** \brief Sort a set of file names in an order that is meaningful for loading them into an mitk::Image. \warning This method assumes that input files are similar in basic properties such as slice thicknes, image orientation, pixel spacing, rows, columns. It should always be ok to put the result of a call to GetSeries(..) into this method. Sorting order is determined by 1. image position along its normal (distance from world origin) 2. acquisition time If P denotes a position and T denotes a time step, this method will order slices from three timesteps like this: \verbatim P1T1 P1T2 P1T3 P2T1 P2T2 P2T3 P3T1 P3T2 P3T3 \endverbatim */ static StringContainer SortSeriesSlices(const StringContainer &unsortedFilenames); public: /** \brief Checks if a specific file is a Philips3D ultrasound DICOM file. */ static bool IsPhilips3DDicom(const std::string &filename); static std::string ReaderImplementationLevelToString( const ReaderImplementationLevel& enumValue ); static std::string PixelSpacingInterpretationToString( const PixelSpacingInterpretation& enumValue ); protected: /** \brief Read a Philips3D ultrasound DICOM file and put into an mitk::Image. */ static bool ReadPhilips3DDicom(const std::string &filename, mitk::Image::Pointer output_image); /** \brief Construct a UID that takes into account sorting criteria from GetSeries(). */ static std::string CreateMoreUniqueSeriesIdentifier( gdcm::Scanner::TagToValue& tagValueMap ); /** \brief Helper for CreateMoreUniqueSeriesIdentifier */ static std::string CreateSeriesIdentifierPart( gdcm::Scanner::TagToValue& tagValueMap, const gdcm::Tag& tag ); /** \brief Helper for CreateMoreUniqueSeriesIdentifier */ static std::string IDifyTagValue(const std::string& value); typedef itk::GDCMImageIO DcmIoType; /** \brief Progress callback for DicomSeriesReader. */ class CallbackCommand : public itk::Command { public: CallbackCommand(UpdateCallBackMethod callback) : m_Callback(callback) { } void Execute(const itk::Object *caller, const itk::EventObject&) { (*this->m_Callback)(static_cast(caller)->GetProgress()); } void Execute(itk::Object *caller, const itk::EventObject&) { (*this->m_Callback)(static_cast(caller)->GetProgress()); } protected: UpdateCallBackMethod m_Callback; }; static void FixSpacingInformation( Image* image, const ImageBlockDescriptor& imageBlockDescriptor ); /** \brief Scan for slice image information */ static void ScanForSliceInformation( const StringContainer &filenames, gdcm::Scanner& scanner ); /** \brief Performs actual loading of a series and creates an image having the specified pixel type. */ - template static void LoadDicom(const StringContainer &filenames, DataNode &node, bool sort, bool check_4d, bool correctTilt, UpdateCallBackMethod callback, Image::Pointer preLoadedImageBlock); /** \brief Feed files into itk::ImageSeriesReader and retrieve a 3D MITK image. \param command can be used for progress reporting */ template static Image::Pointer - LoadDICOMByITK( const StringContainer&, bool correctTilt, const GantryTiltInformation& tiltInfo, DcmIoType::Pointer& io, CallbackCommand* command = NULL, Image::Pointer preLoadedImageBlock = NULL); + LoadDICOMByITK( const StringContainer&, bool correctTilt, const GantryTiltInformation& tiltInfo, DcmIoType::Pointer& io, CallbackCommand* command, Image::Pointer preLoadedImageBlock); + + static + Image::Pointer MultiplexLoadDICOMByITK(const StringContainer&, bool correctTilt, const GantryTiltInformation& tiltInfo, DcmIoType::Pointer& io, CallbackCommand* command, Image::Pointer preLoadedImageBlock); + + static + Image::Pointer MultiplexLoadDICOMByITKScalar(const StringContainer&, bool correctTilt, const GantryTiltInformation& tiltInfo, DcmIoType::Pointer& io, CallbackCommand* command, Image::Pointer preLoadedImageBlock); + + static + Image::Pointer MultiplexLoadDICOMByITKRGBPixel(const StringContainer&, bool correctTilt, const GantryTiltInformation& tiltInfo, DcmIoType::Pointer& io, CallbackCommand* command, Image::Pointer preLoadedImageBlock); + + + template + static + Image::Pointer + LoadDICOMByITK4D( std::list& imageBlocks, ImageBlockDescriptor imageBlockDescriptor, bool correctTilt, const GantryTiltInformation& tiltInfo, DcmIoType::Pointer& io, CallbackCommand* command, Image::Pointer preLoadedImageBlock); + + static + Image::Pointer + MultiplexLoadDICOMByITK4D( std::list& imageBlocks, ImageBlockDescriptor imageBlockDescriptor, bool correctTilt, const GantryTiltInformation& tiltInfo, DcmIoType::Pointer& io, CallbackCommand* command, Image::Pointer preLoadedImageBlock); + + static + Image::Pointer + MultiplexLoadDICOMByITK4DScalar( std::list& imageBlocks, ImageBlockDescriptor imageBlockDescriptor, bool correctTilt, const GantryTiltInformation& tiltInfo, DcmIoType::Pointer& io, CallbackCommand* command, Image::Pointer preLoadedImageBlock); + + static + Image::Pointer + MultiplexLoadDICOMByITK4DRGBPixel( std::list& imageBlocks, ImageBlockDescriptor imageBlockDescriptor, bool correctTilt, const GantryTiltInformation& tiltInfo, DcmIoType::Pointer& io, CallbackCommand* command, Image::Pointer preLoadedImageBlock); + + /** \brief Sort files into time step blocks of a 3D+t image. Called by LoadDicom. Expects to be fed a single list of filenames that have been sorted by GetSeries previously (one map entry). This method will check how many timestep can be filled with given files. Assumption is that the number of time steps is determined by how often the first position in space repeats. I.e. if the first three files in the input parameter all describe the same location in space, we'll construct three lists of files. and sort the remaining files into them. \todo We can probably remove this method if we somehow transfer 3D+t information from GetSeries to LoadDicomSeries. */ static std::list SortIntoBlocksFor3DplusT( const StringContainer& presortedFilenames, const gdcm::Scanner::MappingType& tagValueMappings_, bool sort, bool& canLoadAs4D); /** \brief Defines spatial sorting for sorting by GDCM 2. Sorts by image position along image normal (distance from world origin). In cases of conflict, acquisition time is used as a secondary sort criterium. */ static bool GdcmSortFunction(const gdcm::DataSet &ds1, const gdcm::DataSet &ds2); /** \brief Copy information about files and DICOM tags from ITK's MetaDataDictionary and from the list of input files to the PropertyList of mitk::Image. \todo Tag copy must follow; image level will cause some additional files parsing, probably. */ static void CopyMetaDataToImageProperties( StringContainer filenames, const gdcm::Scanner::MappingType& tagValueMappings_, DcmIoType* io, const ImageBlockDescriptor& blockInfo, Image* image); static void CopyMetaDataToImageProperties( std::list imageBlock, const gdcm::Scanner::MappingType& tagValueMappings_, DcmIoType* io, const ImageBlockDescriptor& blockInfo, Image* image); /** \brief Map between DICOM tags and MITK properties. Uses as a positive list for copying specified DICOM tags (from ITK's ImageIO) to MITK properties. ITK provides MetaDataDictionary entries of form "gggg|eeee" (g = group, e = element), e.g. "0028,0109" (Largest Pixel in Series), which we want to sort as dicom.series.largest_pixel_in_series". */ static const TagToPropertyMapType& GetDICOMTagsToMITKPropertyMap(); }; } #endif /* mitkDicomSeriesReader_h */ diff --git a/Core/Code/IO/mitkDicomSeriesReader.txx b/Core/Code/IO/mitkDicomSeriesReader.txx index e28fa06880..ec4058551c 100644 --- a/Core/Code/IO/mitkDicomSeriesReader.txx +++ b/Core/Code/IO/mitkDicomSeriesReader.txx @@ -1,568 +1,296 @@ /*=================================================================== The Medical Imaging Interaction Toolkit (MITK) Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See LICENSE.txt or http://www.mitk.org for details. ===================================================================*/ #ifndef MITKDICOMSERIESREADER_TXX_ #define MITKDICOMSERIESREADER_TXX_ #include #include #include #include #include #include #include #include namespace mitk { template -void DicomSeriesReader::LoadDicom(const StringContainer &filenames, DataNode &node, bool sort, bool load4D, bool correctTilt, UpdateCallBackMethod callback, Image::Pointer preLoadedImageBlock) +Image::Pointer DicomSeriesReader::LoadDICOMByITK4D( std::list& imageBlocks, ImageBlockDescriptor imageBlockDescriptor, bool correctTilt, const GantryTiltInformation& tiltInfo, DcmIoType::Pointer& io, CallbackCommand* command, Image::Pointer preLoadedImageBlock ) { - const char* previousCLocale = setlocale(LC_NUMERIC, NULL); - setlocale(LC_NUMERIC, "C"); - std::locale previousCppLocale( std::cin.getloc() ); - std::locale l( "C" ); - std::cin.imbue(l); - - ImageBlockDescriptor imageBlockDescriptor; - - const gdcm::Tag tagImagePositionPatient(0x0020,0x0032); // Image Position (Patient) - const gdcm::Tag tagImageOrientation(0x0020, 0x0037); // Image Orientation - const gdcm::Tag tagSeriesInstanceUID(0x0020, 0x000e); // Series Instance UID - const gdcm::Tag tagSOPClassUID(0x0008, 0x0016); // SOP class UID - const gdcm::Tag tagModality(0x0008, 0x0060); // modality - const gdcm::Tag tagPixelSpacing(0x0028, 0x0030); // pixel spacing - const gdcm::Tag tagImagerPixelSpacing(0x0018, 0x1164); // imager pixel spacing - const gdcm::Tag tagNumberOfFrames(0x0028, 0x0008); // number of frames - - try - { - Image::Pointer image = preLoadedImageBlock.IsNull() ? Image::New() : preLoadedImageBlock; - CallbackCommand *command = callback ? new CallbackCommand(callback) : 0; - bool initialize_node = false; - - /* special case for Philips 3D+t ultrasound images */ - if ( DicomSeriesReader::IsPhilips3DDicom(filenames.front().c_str()) ) - { - // TODO what about imageBlockDescriptor? - // TODO what about preLoadedImageBlock? - ReadPhilips3DDicom(filenames.front().c_str(), image); - initialize_node = true; - } - else - { - /* default case: assume "normal" image blocks, possibly 3D+t */ - bool canLoadAs4D(true); - gdcm::Scanner scanner; - ScanForSliceInformation(filenames, scanner); + mitk::Image::Pointer image = mitk::Image::New(); - // need non-const access for map - gdcm::Scanner::MappingType& tagValueMappings = const_cast(scanner.GetMappings()); + // It is 3D+t! Read it and store into mitk image + typedef itk::Image ImageType; + typedef itk::ImageSeriesReader ReaderType; - std::list imageBlocks = SortIntoBlocksFor3DplusT( filenames, tagValueMappings, sort, canLoadAs4D ); - unsigned int volume_count = imageBlocks.size(); + typename ReaderType::Pointer reader = ReaderType::New(); - imageBlockDescriptor.SetSeriesInstanceUID( DicomSeriesReader::ConstCharStarToString( scanner.GetValue( filenames.front().c_str(), tagSeriesInstanceUID ) ) ); - imageBlockDescriptor.SetSOPClassUID( DicomSeriesReader::ConstCharStarToString( scanner.GetValue( filenames.front().c_str(), tagSOPClassUID ) ) ); - imageBlockDescriptor.SetModality( DicomSeriesReader::ConstCharStarToString( scanner.GetValue( filenames.front().c_str(), tagModality ) ) ); - imageBlockDescriptor.SetNumberOfFrames( ConstCharStarToString( scanner.GetValue( filenames.front().c_str(), tagNumberOfFrames ) ) ); - imageBlockDescriptor.SetPixelSpacingInformation( ConstCharStarToString( scanner.GetValue( filenames.front().c_str(), tagPixelSpacing ) ), - ConstCharStarToString( scanner.GetValue( filenames.front().c_str(), tagImagerPixelSpacing ) ) ); + reader->SetImageIO(io); + reader->ReverseOrderOff(); - GantryTiltInformation tiltInfo; + if (command) + { + reader->AddObserver(itk::ProgressEvent(), command); + } - // check possibility of a single slice with many timesteps. In this case, don't check for tilt, no second slice possible - if ( !imageBlocks.empty() && imageBlocks.front().size() > 1 && correctTilt) - { - // check tiltedness here, potentially fixup ITK's loading result by shifting slice contents - // check first and last position slice from tags, make some calculations to detect tilt + unsigned int act_volume = 1u; - std::string firstFilename(imageBlocks.front().front()); - // calculate from first and last slice to minimize rounding errors - std::string secondFilename(imageBlocks.front().back()); + if (preLoadedImageBlock.IsNull()) + { + reader->SetFileNames(imageBlocks.front()); - std::string imagePosition1( ConstCharStarToString( tagValueMappings[ firstFilename.c_str() ][ tagImagePositionPatient ] ) ); - std::string imageOrientation( ConstCharStarToString( tagValueMappings[ firstFilename.c_str() ][ tagImageOrientation ] ) ); - std::string imagePosition2( ConstCharStarToString( tagValueMappings[secondFilename.c_str() ][ tagImagePositionPatient ] ) ); + reader->Update(); - bool ignoredConversionError(-42); // hard to get here, no graceful way to react - Point3D origin1( DICOMStringToPoint3D( imagePosition1, ignoredConversionError ) ); - Point3D origin2( DICOMStringToPoint3D( imagePosition2, ignoredConversionError ) ); + typename ImageType::Pointer readVolume = reader->GetOutput(); + // if we detected that the images are from a tilted gantry acquisition, we need to push some pixels into the right position + if (correctTilt) + { + readVolume = InPlaceFixUpTiltedGeometry( reader->GetOutput(), tiltInfo ); + } - Vector3D right; right.Fill(0.0); - Vector3D up; right.Fill(0.0); // might be down as well, but it is just a name at this point - DICOMStringToOrientationVectors( imageOrientation, right, up, ignoredConversionError ); + unsigned int volume_count = imageBlocks.size(); + image->InitializeByItk( readVolume.GetPointer(), 1, volume_count); + image->SetImportVolume( readVolume->GetBufferPointer(), 0u); - tiltInfo = GantryTiltInformation ( origin1, origin2, right, up, filenames.size()-1 ); - correctTilt = tiltInfo.IsSheared() && tiltInfo.IsRegularGantryTilt(); - } - else - { - correctTilt = false; // we CANNOT do that - } + FixSpacingInformation( image, imageBlockDescriptor ); + } + else + { + StringContainer fakeList; + fakeList.push_back( imageBlocks.front().front() ); + reader->SetFileNames(fakeList); // only ONE first filename to get MetaDataDictionary - imageBlockDescriptor.SetHasGantryTiltCorrected( correctTilt ); + image = preLoadedImageBlock; + } - if (volume_count == 1 || !canLoadAs4D || !load4D) - { + gdcm::Scanner scanner; + ScanForSliceInformation(imageBlockDescriptor.GetFilenames(), scanner); + CopyMetaDataToImageProperties( imageBlocks, scanner.GetMappings(), io, imageBlockDescriptor, image); - DcmIoType::Pointer io; - image = LoadDICOMByITK( imageBlocks.front(), correctTilt, tiltInfo, io, command, preLoadedImageBlock ); // load first 3D block + MITK_DEBUG << "Volume dimension: [" << image->GetDimension(0) << ", " + << image->GetDimension(1) << ", " + << image->GetDimension(2) << ", " + << image->GetDimension(3) << "]"; - imageBlockDescriptor.AddFiles(imageBlocks.front()); // only the first part is loaded - imageBlockDescriptor.SetHasMultipleTimePoints( false ); + MITK_DEBUG << "Volume spacing: [" << image->GetGeometry()->GetSpacing()[0] << ", " + << image->GetGeometry()->GetSpacing()[1] << ", " + << image->GetGeometry()->GetSpacing()[2] << "]"; - FixSpacingInformation( image, imageBlockDescriptor ); - CopyMetaDataToImageProperties( imageBlocks.front(), scanner.GetMappings(), io, imageBlockDescriptor, image); + if ( preLoadedImageBlock.IsNull() ) + { + for (std::list::iterator df_it = ++imageBlocks.begin(); df_it != imageBlocks.end(); ++df_it) + { + reader->SetFileNames(*df_it); + reader->Update(); + typename ImageType::Pointer readVolume = reader->GetOutput(); - initialize_node = true; - } - else if (volume_count > 1) + if (correctTilt) { - imageBlockDescriptor.AddFiles(filenames); // all is loaded - imageBlockDescriptor.SetHasMultipleTimePoints( true ); - // It is 3D+t! Read it and store into mitk image - typedef itk::Image ImageType; - typedef itk::ImageSeriesReader ReaderType; - - DcmIoType::Pointer io = DcmIoType::New(); - typename ReaderType::Pointer reader = ReaderType::New(); - - reader->SetImageIO(io); - reader->ReverseOrderOff(); - - if (command) - { - reader->AddObserver(itk::ProgressEvent(), command); - } - - unsigned int act_volume = 1u; - - if (preLoadedImageBlock.IsNull()) - { - reader->SetFileNames(imageBlocks.front()); - - reader->Update(); - - typename ImageType::Pointer readVolume = reader->GetOutput(); - // if we detected that the images are from a tilted gantry acquisition, we need to push some pixels into the right position - if (correctTilt) - { - readVolume = InPlaceFixUpTiltedGeometry( reader->GetOutput(), tiltInfo ); - } - - image->InitializeByItk( readVolume.GetPointer(), 1, volume_count); - image->SetImportVolume( readVolume->GetBufferPointer(), 0u); - - FixSpacingInformation( image, imageBlockDescriptor ); - } - else - { - StringContainer fakeList; - fakeList.push_back( imageBlocks.front().front() ); - reader->SetFileNames(fakeList); // only ONE first filename to get MetaDataDictionary - } - - gdcm::Scanner scanner; - ScanForSliceInformation(filenames, scanner); - CopyMetaDataToImageProperties( imageBlocks, scanner.GetMappings(), io, imageBlockDescriptor, image); - - MITK_DEBUG << "Volume dimension: [" << image->GetDimension(0) << ", " - << image->GetDimension(1) << ", " - << image->GetDimension(2) << ", " - << image->GetDimension(3) << "]"; - - MITK_DEBUG << "Volume spacing: [" << image->GetGeometry()->GetSpacing()[0] << ", " - << image->GetGeometry()->GetSpacing()[1] << ", " - << image->GetGeometry()->GetSpacing()[2] << "]"; - - if ( preLoadedImageBlock.IsNull() ) - { - for (std::list::iterator df_it = ++imageBlocks.begin(); df_it != imageBlocks.end(); ++df_it) - { - reader->SetFileNames(*df_it); - reader->Update(); - typename ImageType::Pointer readVolume = reader->GetOutput(); - - if (correctTilt) - { - readVolume = InPlaceFixUpTiltedGeometry( reader->GetOutput(), tiltInfo ); - } - - image->SetImportVolume(readVolume->GetBufferPointer(), act_volume++); - } - } - - initialize_node = true; + readVolume = InPlaceFixUpTiltedGeometry( reader->GetOutput(), tiltInfo ); } + image->SetImportVolume(readVolume->GetBufferPointer(), act_volume++); } - - if (initialize_node) - { - // forward some image properties to node - node.GetPropertyList()->ConcatenatePropertyList( image->GetPropertyList(), true ); - - node.SetData( image ); - setlocale(LC_NUMERIC, previousCLocale); - std::cin.imbue(previousCppLocale); - } - - MITK_DEBUG << "--------------------------------------------------------------------------------"; - MITK_DEBUG << "DICOM files loaded (from series UID " << imageBlockDescriptor.GetSeriesInstanceUID() << "):"; - MITK_DEBUG << " " << imageBlockDescriptor.GetFilenames().size() << " '" << imageBlockDescriptor.GetModality() << "' files (" << imageBlockDescriptor.GetSOPClassUIDAsString() << ") loaded into 1 mitk::Image"; - MITK_DEBUG << " multi-frame: " << (imageBlockDescriptor.IsMultiFrameImage()?"Yes":"No"); - MITK_DEBUG << " reader support: " << ReaderImplementationLevelToString(imageBlockDescriptor.GetReaderImplementationLevel()); - MITK_DEBUG << " pixel spacing type: " << PixelSpacingInterpretationToString( imageBlockDescriptor.GetPixelSpacingType() ) << " " << image->GetGeometry()->GetSpacing()[0] << "/" << image->GetGeometry()->GetSpacing()[0]; - MITK_DEBUG << " gantry tilt corrected: " << (imageBlockDescriptor.HasGantryTiltCorrected()?"Yes":"No"); - MITK_DEBUG << " 3D+t: " << (imageBlockDescriptor.HasMultipleTimePoints()?"Yes":"No"); - MITK_DEBUG << "--------------------------------------------------------------------------------"; } - catch (std::exception& e) - { - // reset locale then throw up - setlocale(LC_NUMERIC, previousCLocale); - std::cin.imbue(previousCppLocale); - - MITK_DEBUG << "Caught exception in DicomSeriesReader::LoadDicom"; - throw e; - } + return image; } template Image::Pointer DicomSeriesReader::LoadDICOMByITK( const StringContainer& filenames, bool correctTilt, const GantryTiltInformation& tiltInfo, DcmIoType::Pointer& io, CallbackCommand* command, Image::Pointer preLoadedImageBlock ) { /******** Normal Case, 3D (also for GDCM < 2 usable) ***************/ mitk::Image::Pointer image = mitk::Image::New(); typedef itk::Image ImageType; typedef itk::ImageSeriesReader ReaderType; io = DcmIoType::New(); typename ReaderType::Pointer reader = ReaderType::New(); reader->SetImageIO(io); reader->ReverseOrderOff(); if (command) { reader->AddObserver(itk::ProgressEvent(), command); } if (preLoadedImageBlock.IsNull()) { reader->SetFileNames(filenames); reader->Update(); typename ImageType::Pointer readVolume = reader->GetOutput(); // if we detected that the images are from a tilted gantry acquisition, we need to push some pixels into the right position if (correctTilt) { readVolume = InPlaceFixUpTiltedGeometry( reader->GetOutput(), tiltInfo ); } image->InitializeByItk(readVolume.GetPointer()); image->SetImportVolume(readVolume->GetBufferPointer()); } else { image = preLoadedImageBlock; StringContainer fakeList; fakeList.push_back( filenames.front() ); reader->SetFileNames( fakeList ); // we always need to load at least one file to get the MetaDataDictionary reader->Update(); } MITK_DEBUG << "Volume dimension: [" << image->GetDimension(0) << ", " << image->GetDimension(1) << ", " << image->GetDimension(2) << "]"; MITK_DEBUG << "Volume spacing: [" << image->GetGeometry()->GetSpacing()[0] << ", " << image->GetGeometry()->GetSpacing()[1] << ", " << image->GetGeometry()->GetSpacing()[2] << "]"; return image; } -void -DicomSeriesReader::ScanForSliceInformation(const StringContainer &filenames, gdcm::Scanner& scanner) -{ - const gdcm::Tag tagImagePositionPatient(0x0020,0x0032); //Image position (Patient) - scanner.AddTag(tagImagePositionPatient); - - const gdcm::Tag tagSeriesInstanceUID(0x0020, 0x000e); // Series Instance UID - scanner.AddTag(tagSeriesInstanceUID); - - const gdcm::Tag tagImageOrientation(0x0020,0x0037); //Image orientation - scanner.AddTag(tagImageOrientation); - - const gdcm::Tag tagSliceLocation(0x0020, 0x1041); // slice location - scanner.AddTag( tagSliceLocation ); - - const gdcm::Tag tagInstanceNumber(0x0020, 0x0013); // (image) instance number - scanner.AddTag( tagInstanceNumber ); - - const gdcm::Tag tagSOPInstanceNumber(0x0008, 0x0018); // SOP instance number - scanner.AddTag( tagSOPInstanceNumber ); - - const gdcm::Tag tagPixelSpacing(0x0028, 0x0030); // Pixel Spacing - scanner.AddTag( tagPixelSpacing ); - - const gdcm::Tag tagImagerPixelSpacing(0x0018, 0x1164); // Imager Pixel Spacing - scanner.AddTag( tagImagerPixelSpacing ); - - const gdcm::Tag tagModality(0x0008, 0x0060); // Modality - scanner.AddTag( tagModality ); - - const gdcm::Tag tagSOPClassUID(0x0008, 0x0016); // SOP Class UID - scanner.AddTag( tagSOPClassUID ); - - const gdcm::Tag tagNumberOfFrames(0x0028, 0x0008); // number of frames - scanner.AddTag( tagNumberOfFrames ); - - scanner.Scan(filenames); // make available image information for each file -} - -std::list -DicomSeriesReader::SortIntoBlocksFor3DplusT( - const StringContainer& presortedFilenames, - const gdcm::Scanner::MappingType& tagValueMappings, - bool /*sort*/, - bool& canLoadAs4D ) -{ - std::list imageBlocks; - - // ignore sort request, because most likely re-sorting is now needed due to changes in GetSeries(bug #8022) - StringContainer sorted_filenames = DicomSeriesReader::SortSeriesSlices(presortedFilenames); - - std::string firstPosition; - unsigned int numberOfBlocks(0); // number of 3D image blocks - - static const gdcm::Tag tagImagePositionPatient(0x0020,0x0032); //Image position (Patient) - - // loop files to determine number of image blocks - for (StringContainer::const_iterator fileIter = sorted_filenames.begin(); - fileIter != sorted_filenames.end(); - ++fileIter) - { - gdcm::Scanner::TagToValue tagToValueMap = tagValueMappings.find( fileIter->c_str() )->second; - - if(tagToValueMap.find(tagImagePositionPatient) == tagToValueMap.end()) - { - // we expect to get images w/ missing position information ONLY as separated blocks. - assert( presortedFilenames.size() == 1 ); - numberOfBlocks = 1; - break; - } - - std::string position = tagToValueMap.find(tagImagePositionPatient)->second; - MITK_DEBUG << " " << *fileIter << " at " << position; - if (firstPosition.empty()) - { - firstPosition = position; - } - - if ( position == firstPosition ) - { - ++numberOfBlocks; - } - else - { - break; // enough information to know the number of image blocks - } - } - - MITK_DEBUG << " ==> Assuming " << numberOfBlocks << " time steps"; - - if (numberOfBlocks == 0) return imageBlocks; // only possible if called with no files - - - // loop files to sort them into image blocks - unsigned int numberOfExpectedSlices(0); - for (unsigned int block = 0; block < numberOfBlocks; ++block) - { - StringContainer filesOfCurrentBlock; - - for ( StringContainer::const_iterator fileIter = sorted_filenames.begin() + block; - fileIter != sorted_filenames.end(); - //fileIter += numberOfBlocks) // TODO shouldn't this work? give invalid iterators on first attempts - ) - { - filesOfCurrentBlock.push_back( *fileIter ); - for (unsigned int b = 0; b < numberOfBlocks; ++b) - { - if (fileIter != sorted_filenames.end()) - ++fileIter; - } - } - - imageBlocks.push_back(filesOfCurrentBlock); - - if (block == 0) - { - numberOfExpectedSlices = filesOfCurrentBlock.size(); - } - else - { - if (filesOfCurrentBlock.size() != numberOfExpectedSlices) - { - MITK_WARN << "DicomSeriesReader expected " << numberOfBlocks - << " image blocks of " - << numberOfExpectedSlices - << " images each. Block " - << block - << " got " - << filesOfCurrentBlock.size() - << " instead. Cannot load this as 3D+t"; // TODO implement recovery (load as many slices 3D+t as much as possible) - canLoadAs4D = false; - } - } - } - - return imageBlocks; -} - template typename ImageType::Pointer DicomSeriesReader::InPlaceFixUpTiltedGeometry( ImageType* input, const GantryTiltInformation& tiltInfo ) { typedef itk::ResampleImageFilter ResampleFilterType; typename ResampleFilterType::Pointer resampler = ResampleFilterType::New(); resampler->SetInput( input ); /* Transform for a point is - transform from actual position to index coordinates - apply a shear that undoes the gantry tilt - transform back into world coordinates Anybody who does this in a simpler way: don't forget to write up how and why your solution works */ typedef itk::ScalableAffineTransform< double, ImageType::ImageDimension > TransformType; typename TransformType::Pointer transformShear = TransformType::New(); /** - apply a shear and spacing correction to the image block that corrects the ITK reader's error - ITK ignores the shear and loads slices into an orthogonal volume - ITK calculates the spacing from the origin distance, which is more than the actual spacing with gantry tilt images - to undo the effect - we have calculated some information in tiltInfo: - the shift in Y direction that is added with each additional slice is the most important information - the Y-shift is calculated in mm world coordinates - we apply a shearing transformation to the ITK-read image volume - to do this locally, - we transform the image volume back to origin and "normal" orientation by applying the inverse of its transform (this brings us into the image's "index coordinate" system) - we apply a shear with the Y-shift factor put into a unit transform at row 1, col 2 - we transform the image volume back to its actual position (from index to world coordinates) - we lastly apply modify the image spacing in z direction by replacing this number with the correctly calulcated inter-slice distance */ ScalarType factor = tiltInfo.GetMatrixCoefficientForCorrectionInWorldCoordinates() / input->GetSpacing()[1]; // row 1, column 2 corrects shear in parallel to Y axis, proportional to distance in Z direction transformShear->Shear( 1, 2, factor ); typename TransformType::Pointer imageIndexToWorld = TransformType::New(); imageIndexToWorld->SetOffset( input->GetOrigin().GetVectorFromOrigin() ); typename TransformType::MatrixType indexToWorldMatrix; indexToWorldMatrix = input->GetDirection(); typename ImageType::DirectionType scale; for ( unsigned int i = 0; i < ImageType::ImageDimension; i++ ) { scale[i][i] = input->GetSpacing()[i]; } indexToWorldMatrix *= scale; imageIndexToWorld->SetMatrix( indexToWorldMatrix ); typename TransformType::Pointer imageWorldToIndex = TransformType::New(); imageIndexToWorld->GetInverse( imageWorldToIndex ); typename TransformType::Pointer gantryTiltCorrection = TransformType::New(); gantryTiltCorrection->Compose( imageWorldToIndex ); gantryTiltCorrection->Compose( transformShear ); gantryTiltCorrection->Compose( imageIndexToWorld ); resampler->SetTransform( gantryTiltCorrection ); typedef itk::LinearInterpolateImageFunction< ImageType, double > InterpolatorType; typename InterpolatorType::Pointer interpolator = InterpolatorType::New(); resampler->SetInterpolator( interpolator ); /* This would be the right place to invent a meaningful value for positions outside of the image. For CT, HU -1000 might be meaningful, but a general solution seems not possible. Even for CT, -1000 would only look natural for many not all images. */ // TODO use (0028,0120) Pixel Padding Value if present resampler->SetDefaultPixelValue( itk::NumericTraits< typename ImageType::PixelType >::min() ); // adjust size in Y direction! (maybe just transform the outer last pixel to see how much space we would need resampler->SetOutputParametersFromImage( input ); // we basically need the same image again, just sheared // if tilt positive, then we need additional pixels BELOW origin, otherwise we need pixels behind the end of the block // in any case we need more size to accomodate shifted slices typename ImageType::SizeType largerSize = resampler->GetSize(); // now the resampler already holds the input image's size. largerSize[1] += static_cast(tiltInfo.GetTiltCorrectedAdditionalSize() / input->GetSpacing()[1]+ 2.0); resampler->SetSize( largerSize ); // in SOME cases this additional size is below/behind origin if ( tiltInfo.GetMatrixCoefficientForCorrectionInWorldCoordinates() > 0.0 ) { typename ImageType::DirectionType imageDirection = input->GetDirection(); Vector3D yDirection; yDirection[0] = imageDirection[0][1]; yDirection[1] = imageDirection[1][1]; yDirection[2] = imageDirection[2][1]; yDirection.Normalize(); typename ImageType::PointType shiftedOrigin; shiftedOrigin = input->GetOrigin(); // add some pixels to make everything fit shiftedOrigin[0] -= yDirection[0] * (tiltInfo.GetTiltCorrectedAdditionalSize() + 1.0 * input->GetSpacing()[1]); shiftedOrigin[1] -= yDirection[1] * (tiltInfo.GetTiltCorrectedAdditionalSize() + 1.0 * input->GetSpacing()[1]); shiftedOrigin[2] -= yDirection[2] * (tiltInfo.GetTiltCorrectedAdditionalSize() + 1.0 * input->GetSpacing()[1]); resampler->SetOutputOrigin( shiftedOrigin ); } resampler->Update(); typename ImageType::Pointer result = resampler->GetOutput(); // ImageSeriesReader calculates z spacing as the distance between the first two origins. // This is not correct in case of gantry tilt, so we set our calculated spacing. typename ImageType::SpacingType correctedSpacing = result->GetSpacing(); correctedSpacing[2] = tiltInfo.GetRealZSpacing(); result->SetSpacing( correctedSpacing ); return result; } - } #endif diff --git a/Core/Code/files.cmake b/Core/Code/files.cmake index 73f6320348..c8a588c37b 100644 --- a/Core/Code/files.cmake +++ b/Core/Code/files.cmake @@ -1,399 +1,407 @@ set(H_FILES Algorithms/itkImportMitkImageContainer.h Algorithms/itkImportMitkImageContainer.txx Algorithms/itkLocalVariationImageFilter.h Algorithms/itkLocalVariationImageFilter.txx Algorithms/itkMITKScalarImageToHistogramGenerator.h Algorithms/itkMITKScalarImageToHistogramGenerator.txx Algorithms/itkTotalVariationDenoisingImageFilter.h Algorithms/itkTotalVariationDenoisingImageFilter.txx Algorithms/itkTotalVariationSingleIterationImageFilter.h Algorithms/itkTotalVariationSingleIterationImageFilter.txx Algorithms/mitkBilateralFilter.h Algorithms/mitkBilateralFilter.cpp Algorithms/mitkInstantiateAccessFunctions.h Algorithms/mitkPixelTypeList.h Algorithms/mitkPPArithmeticDec.h Algorithms/mitkPPArgCount.h Algorithms/mitkPPCat.h Algorithms/mitkPPConfig.h Algorithms/mitkPPControlExprIIf.h Algorithms/mitkPPControlIf.h Algorithms/mitkPPControlIIf.h Algorithms/mitkPPDebugError.h Algorithms/mitkPPDetailAutoRec.h Algorithms/mitkPPDetailDMCAutoRec.h Algorithms/mitkPPExpand.h Algorithms/mitkPPFacilitiesEmpty.h Algorithms/mitkPPFacilitiesExpand.h Algorithms/mitkPPLogicalBool.h Algorithms/mitkPPRepetitionDetailDMCFor.h Algorithms/mitkPPRepetitionDetailEDGFor.h Algorithms/mitkPPRepetitionDetailFor.h Algorithms/mitkPPRepetitionDetailMSVCFor.h Algorithms/mitkPPRepetitionFor.h Algorithms/mitkPPSeqElem.h Algorithms/mitkPPSeqForEach.h Algorithms/mitkPPSeqForEachProduct.h Algorithms/mitkPPSeq.h Algorithms/mitkPPSeqEnum.h Algorithms/mitkPPSeqSize.h Algorithms/mitkPPSeqToTuple.h Algorithms/mitkPPStringize.h Algorithms/mitkPPTupleEat.h Algorithms/mitkPPTupleElem.h Algorithms/mitkPPTupleRem.h Algorithms/mitkClippedSurfaceBoundsCalculator.h Algorithms/mitkExtractSliceFilter.h Algorithms/mitkConvert2Dto3DImageFilter.h Algorithms/mitkPlaneClipping.h Common/mitkExceptionMacro.h Common/mitkServiceBaseObject.h Common/mitkTestingMacros.h Common/mitkTesting.h DataManagement/mitkProportionalTimeGeometry.h DataManagement/mitkTimeGeometry.h DataManagement/mitkImageAccessByItk.h DataManagement/mitkImageCast.h DataManagement/mitkImagePixelAccessor.h DataManagement/mitkImagePixelReadAccessor.h DataManagement/mitkImagePixelWriteAccessor.h DataManagement/mitkImageReadAccessor.h DataManagement/mitkImageWriteAccessor.h DataManagement/mitkITKImageImport.h DataManagement/mitkITKImageImport.txx DataManagement/mitkImageToItk.h DataManagement/mitkImageToItk.txx DataManagement/mitkTimeSlicedGeometry.h # Deprecated, empty for compatibilty reasons. Interactions/mitkEventMapperAddOn.h Interfaces/mitkIDataNodeReader.h Rendering/mitkLocalStorageHandler.h IO/mitkPixelTypeTraits.h ) set(CPP_FILES Algorithms/mitkBaseDataSource.cpp Algorithms/mitkCompareImageFilter.cpp Algorithms/mitkDataNodeSource.cpp Algorithms/mitkGeometry2DDataToSurfaceFilter.cpp Algorithms/mitkHistogramGenerator.cpp Algorithms/mitkImageChannelSelector.cpp Algorithms/mitkImageSliceSelector.cpp Algorithms/mitkImageSource.cpp Algorithms/mitkImageTimeSelector.cpp Algorithms/mitkImageToImageFilter.cpp Algorithms/mitkImageToSurfaceFilter.cpp Algorithms/mitkPointSetSource.cpp Algorithms/mitkPointSetToPointSetFilter.cpp Algorithms/mitkRGBToRGBACastImageFilter.cpp Algorithms/mitkSubImageSelector.cpp Algorithms/mitkSurfaceSource.cpp Algorithms/mitkSurfaceToImageFilter.cpp Algorithms/mitkSurfaceToSurfaceFilter.cpp Algorithms/mitkUIDGenerator.cpp Algorithms/mitkVolumeCalculator.cpp Algorithms/mitkClippedSurfaceBoundsCalculator.cpp Algorithms/mitkExtractSliceFilter.cpp Algorithms/mitkConvert2Dto3DImageFilter.cpp Controllers/mitkBaseController.cpp Controllers/mitkCallbackFromGUIThread.cpp Controllers/mitkCameraController.cpp Controllers/mitkCameraRotationController.cpp Controllers/mitkCoreActivator.cpp Controllers/mitkFocusManager.cpp Controllers/mitkLimitedLinearUndo.cpp Controllers/mitkOperationEvent.cpp Controllers/mitkPlanePositionManager.cpp Controllers/mitkProgressBar.cpp Controllers/mitkRenderingManager.cpp Controllers/mitkSliceNavigationController.cpp Controllers/mitkSlicesCoordinator.cpp Controllers/mitkSlicesRotator.cpp Controllers/mitkSlicesSwiveller.cpp Controllers/mitkStatusBar.cpp Controllers/mitkStepper.cpp Controllers/mitkTestManager.cpp Controllers/mitkUndoController.cpp Controllers/mitkVerboseLimitedLinearUndo.cpp Controllers/mitkVtkInteractorCameraController.cpp Controllers/mitkVtkLayerController.cpp DataManagement/mitkProportionalTimeGeometry.cpp DataManagement/mitkTimeGeometry.cpp DataManagement/mitkAbstractTransformGeometry.cpp DataManagement/mitkAnnotationProperty.cpp DataManagement/mitkApplicationCursor.cpp DataManagement/mitkBaseData.cpp DataManagement/mitkBaseProperty.cpp DataManagement/mitkClippingProperty.cpp DataManagement/mitkChannelDescriptor.cpp DataManagement/mitkColorProperty.cpp DataManagement/mitkDataStorage.cpp # DataManagement/mitkDataTree.cpp DataManagement/mitkDataNode.cpp DataManagement/mitkDataNodeFactory.cpp # DataManagement/mitkDataTreeStorage.cpp DataManagement/mitkDisplayGeometry.cpp DataManagement/mitkEnumerationProperty.cpp DataManagement/mitkGeometry2D.cpp DataManagement/mitkGeometry2DData.cpp DataManagement/mitkGeometry3D.cpp DataManagement/mitkGeometryData.cpp DataManagement/mitkGroupTagProperty.cpp DataManagement/mitkImage.cpp DataManagement/mitkImageAccessorBase.cpp DataManagement/mitkImageCaster.cpp DataManagement/mitkImageCastPart1.cpp DataManagement/mitkImageCastPart2.cpp DataManagement/mitkImageCastPart3.cpp DataManagement/mitkImageCastPart4.cpp DataManagement/mitkImageDataItem.cpp DataManagement/mitkImageDescriptor.cpp DataManagement/mitkImageVtkAccessor.cpp DataManagement/mitkImageStatisticsHolder.cpp DataManagement/mitkLandmarkBasedCurvedGeometry.cpp DataManagement/mitkLandmarkProjectorBasedCurvedGeometry.cpp DataManagement/mitkLandmarkProjector.cpp DataManagement/mitkLevelWindow.cpp DataManagement/mitkLevelWindowManager.cpp DataManagement/mitkLevelWindowPreset.cpp DataManagement/mitkLevelWindowProperty.cpp DataManagement/mitkLookupTable.cpp DataManagement/mitkLookupTables.cpp # specializations of GenericLookupTable DataManagement/mitkMemoryUtilities.cpp DataManagement/mitkModalityProperty.cpp DataManagement/mitkModeOperation.cpp DataManagement/mitkNodePredicateAnd.cpp DataManagement/mitkNodePredicateBase.cpp DataManagement/mitkNodePredicateCompositeBase.cpp DataManagement/mitkNodePredicateData.cpp DataManagement/mitkNodePredicateDataType.cpp DataManagement/mitkNodePredicateDimension.cpp DataManagement/mitkNodePredicateFirstLevel.cpp DataManagement/mitkNodePredicateNot.cpp DataManagement/mitkNodePredicateOr.cpp DataManagement/mitkNodePredicateProperty.cpp DataManagement/mitkNodePredicateSource.cpp DataManagement/mitkPlaneOrientationProperty.cpp DataManagement/mitkPlaneGeometry.cpp DataManagement/mitkPlaneOperation.cpp DataManagement/mitkPointOperation.cpp DataManagement/mitkPointSet.cpp DataManagement/mitkProperties.cpp DataManagement/mitkPropertyList.cpp DataManagement/mitkRestorePlanePositionOperation.cpp DataManagement/mitkRotationOperation.cpp DataManagement/mitkSlicedData.cpp DataManagement/mitkSlicedGeometry3D.cpp DataManagement/mitkSmartPointerProperty.cpp DataManagement/mitkStandaloneDataStorage.cpp DataManagement/mitkStateTransitionOperation.cpp DataManagement/mitkStringProperty.cpp DataManagement/mitkSurface.cpp DataManagement/mitkSurfaceOperation.cpp DataManagement/mitkThinPlateSplineCurvedGeometry.cpp DataManagement/mitkTransferFunction.cpp DataManagement/mitkTransferFunctionProperty.cpp DataManagement/mitkTransferFunctionInitializer.cpp DataManagement/mitkVector.cpp DataManagement/mitkVtkInterpolationProperty.cpp DataManagement/mitkVtkRepresentationProperty.cpp DataManagement/mitkVtkResliceInterpolationProperty.cpp DataManagement/mitkVtkScalarModeProperty.cpp DataManagement/mitkVtkVolumeRenderingProperty.cpp DataManagement/mitkWeakPointerProperty.cpp DataManagement/mitkRenderingModeProperty.cpp DataManagement/mitkShaderProperty.cpp DataManagement/mitkResliceMethodProperty.cpp DataManagement/mitkMaterial.cpp DataManagement/mitkPointSetShapeProperty.cpp DataManagement/mitkFloatPropertyExtension.cpp DataManagement/mitkIntPropertyExtension.cpp DataManagement/mitkPropertyExtension.cpp DataManagement/mitkPropertyFilter.cpp DataManagement/mitkPropertyAliases.cpp DataManagement/mitkPropertyDescriptions.cpp DataManagement/mitkPropertyExtensions.cpp DataManagement/mitkPropertyFilters.cpp Interactions/mitkAction.cpp Interactions/mitkAffineInteractor.cpp Interactions/mitkBindDispatcherInteractor.cpp Interactions/mitkCoordinateSupplier.cpp Interactions/mitkDataInteractor.cpp Interactions/mitkDispatcher.cpp Interactions/mitkDisplayCoordinateOperation.cpp Interactions/mitkDisplayInteractor.cpp Interactions/mitkDisplayPositionEvent.cpp # Interactions/mitkDisplayVectorInteractorLevelWindow.cpp # legacy, prob even now unneeded # Interactions/mitkDisplayVectorInteractorScroll.cpp Interactions/mitkEvent.cpp Interactions/mitkEventConfig.cpp Interactions/mitkEventDescription.cpp Interactions/mitkEventFactory.cpp Interactions/mitkInteractionEventHandler.cpp Interactions/mitkEventMapper.cpp Interactions/mitkEventStateMachine.cpp Interactions/mitkGlobalInteraction.cpp Interactions/mitkInteractor.cpp Interactions/mitkInternalEvent.cpp Interactions/mitkInteractionEvent.cpp Interactions/mitkInteractionEventConst.cpp Interactions/mitkInteractionPositionEvent.cpp Interactions/mitkInteractionKeyEvent.cpp Interactions/mitkMousePressEvent.cpp Interactions/mitkMouseMoveEvent.cpp Interactions/mitkMouseReleaseEvent.cpp Interactions/mitkMouseWheelEvent.cpp Interactions/mitkMouseDoubleClickEvent.cpp Interactions/mitkMouseModeSwitcher.cpp Interactions/mitkMouseMovePointSetInteractor.cpp Interactions/mitkMoveBaseDataInteractor.cpp Interactions/mitkNodeDepententPointSetInteractor.cpp Interactions/mitkPointSetDataInteractor.cpp Interactions/mitkPointSetInteractor.cpp Interactions/mitkPositionEvent.cpp Interactions/mitkPositionTracker.cpp Interactions/mitkStateMachineAction.cpp Interactions/mitkStateMachineCondition.cpp Interactions/mitkStateMachineState.cpp Interactions/mitkStateMachineTransition.cpp Interactions/mitkState.cpp Interactions/mitkStateMachineContainer.cpp Interactions/mitkStateEvent.cpp Interactions/mitkStateMachine.cpp Interactions/mitkStateMachineFactory.cpp Interactions/mitkTransition.cpp Interactions/mitkWheelEvent.cpp Interactions/mitkKeyEvent.cpp Interactions/mitkVtkEventAdapter.cpp Interactions/mitkVtkInteractorStyle.cxx Interactions/mitkCrosshairPositionEvent.cpp Interfaces/mitkInteractionEventObserver.cpp Interfaces/mitkIShaderRepository.cpp Interfaces/mitkIPropertyAliases.cpp Interfaces/mitkIPropertyDescriptions.cpp Interfaces/mitkIPropertyExtensions.cpp Interfaces/mitkIPropertyFilters.cpp IO/mitkBaseDataIOFactory.cpp IO/mitkCoreDataNodeReader.cpp IO/mitkDicomSeriesReader.cpp + IO/mitkDicomSR_LoadDICOMScalar.cpp + IO/mitkDicomSR_LoadDICOMScalar4D.cpp + IO/mitkDicomSR_LoadDICOMRGBPixel.cpp + IO/mitkDicomSR_LoadDICOMRGBPixel4D.cpp + IO/mitkDicomSR_ImageBlockDescriptor.cpp + IO/mitkDicomSR_GantryTiltInformation.cpp + IO/mitkDicomSR_SliceGroupingResult.cpp + IO/mitkFileReader.cpp IO/mitkFileSeriesReader.cpp IO/mitkFileWriter.cpp # IO/mitkIpPicGet.c IO/mitkImageGenerator.cpp IO/mitkImageWriter.cpp IO/mitkImageWriterFactory.cpp IO/mitkItkImageFileIOFactory.cpp IO/mitkItkImageFileReader.cpp IO/mitkItkLoggingAdapter.cpp IO/mitkItkPictureWrite.cpp IO/mitkIOUtil.cpp IO/mitkLookupTableProperty.cpp IO/mitkOperation.cpp # IO/mitkPicFileIOFactory.cpp # IO/mitkPicFileReader.cpp # IO/mitkPicFileWriter.cpp # IO/mitkPicHelper.cpp # IO/mitkPicVolumeTimeSeriesIOFactory.cpp # IO/mitkPicVolumeTimeSeriesReader.cpp IO/mitkPixelType.cpp IO/mitkPointSetIOFactory.cpp IO/mitkPointSetReader.cpp IO/mitkPointSetWriter.cpp IO/mitkPointSetWriterFactory.cpp IO/mitkRawImageFileReader.cpp IO/mitkStandardFileLocations.cpp IO/mitkSTLFileIOFactory.cpp IO/mitkSTLFileReader.cpp IO/mitkSurfaceVtkWriter.cpp IO/mitkSurfaceVtkWriterFactory.cpp IO/mitkVtkLoggingAdapter.cpp IO/mitkVtiFileIOFactory.cpp IO/mitkVtiFileReader.cpp IO/mitkVtkImageIOFactory.cpp IO/mitkVtkImageReader.cpp IO/mitkVtkSurfaceIOFactory.cpp IO/mitkVtkSurfaceReader.cpp IO/vtkPointSetXMLParser.cpp IO/mitkLog.cpp Rendering/mitkBaseRenderer.cpp Rendering/mitkVtkMapper.cpp Rendering/mitkRenderWindowFrame.cpp Rendering/mitkGeometry2DDataMapper2D.cpp Rendering/mitkGeometry2DDataVtkMapper3D.cpp Rendering/mitkGLMapper.cpp Rendering/mitkGradientBackground.cpp Rendering/mitkManufacturerLogo.cpp Rendering/mitkMapper.cpp Rendering/mitkPointSetGLMapper2D.cpp Rendering/mitkPointSetVtkMapper2D.cpp Rendering/mitkPointSetVtkMapper3D.cpp Rendering/mitkPolyDataGLMapper2D.cpp Rendering/mitkSurfaceGLMapper2D.cpp Rendering/mitkSurfaceVtkMapper3D.cpp Rendering/mitkVolumeDataVtkMapper3D.cpp Rendering/mitkVtkPropRenderer.cpp Rendering/mitkVtkWidgetRendering.cpp Rendering/vtkMitkRectangleProp.cpp Rendering/vtkMitkRenderProp.cpp Rendering/mitkVtkEventProvider.cpp Rendering/mitkRenderWindow.cpp Rendering/mitkRenderWindowBase.cpp Rendering/mitkShaderRepository.cpp Rendering/mitkImageVtkMapper2D.cpp Rendering/vtkMitkThickSlicesFilter.cpp Rendering/vtkMitkLevelWindowFilter.cpp Rendering/vtkNeverTranslucentTexture.cpp Rendering/mitkRenderingTestHelper.cpp Rendering/mitkOverlay.cpp Rendering/mitkVtkOverlay.cpp Rendering/mitkVtkOverlay2D.cpp Rendering/mitkVtkOverlay3D.cpp Rendering/mitkOverlayManager.cpp Rendering/mitkAbstractOverlayLayouter.cpp Rendering/mitkTextOverlay2D.cpp Rendering/mitkTextOverlay3D.cpp Rendering/mitkLabelOverlay3D.cpp Rendering/mitkOverlay2DLayouter.cpp Common/mitkException.cpp Common/mitkCommon.h Common/mitkCoreObjectFactoryBase.cpp Common/mitkCoreObjectFactory.cpp Common/mitkCoreServices.cpp ) list(APPEND CPP_FILES ${CppMicroServices_SOURCES}) set(RESOURCE_FILES Interactions/globalConfig.xml Interactions/DisplayInteraction.xml Interactions/DisplayConfig.xml Interactions/DisplayConfigPACS.xml Interactions/DisplayConfigPACSPan.xml Interactions/DisplayConfigPACSScroll.xml Interactions/DisplayConfigPACSZoom.xml Interactions/DisplayConfigPACSLevelWindow.xml Interactions/DisplayConfigMITK.xml Interactions/PointSet.xml Interactions/Legacy/StateMachine.xml Interactions/Legacy/DisplayConfigMITKTools.xml Interactions/PointSetConfig.xml Shaders/mitkShaderLighting.xml mitkLevelWindowPresets.xml )