diff --git a/Core/Code/IO/mitkDicomSeriesReader.cpp b/Core/Code/IO/mitkDicomSeriesReader.cpp index d2ef6c7c26..543760d7d7 100644 --- a/Core/Code/IO/mitkDicomSeriesReader.cpp +++ b/Core/Code/IO/mitkDicomSeriesReader.cpp @@ -1,1327 +1,1323 @@ /*=================================================================== 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 "mitkProperties.h" namespace mitk { typedef itk::GDCMSeriesFileNames DcmFileNamesGeneratorType; 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; } -Vector3D DicomSeriesReader::SliceGroupingAnalysisResult::GetInterSliceOffset() -{ - return m_InterSliceOffset; -} - 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::FlagGantryTilt(Vector3D interSliceOffset) +void DicomSeriesReader::SliceGroupingAnalysisResult::FlagGantryTilt() { - m_InterSliceOffset = interSliceOffset; 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"; initialized = true; } return dictionary; } DataNode::Pointer DicomSeriesReader::LoadDicomSeries(const StringContainer &filenames, bool sort, bool check_4d, bool correctTilt, UpdateCallBackMethod callback) { DataNode::Pointer node = DataNode::New(); if (DicomSeriesReader::LoadDicomSeries(filenames, *node, sort, check_4d, correctTilt, callback)) { 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) { if( filenames.empty() ) { MITK_WARN << "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(); switch (io->GetComponentType()) { case DcmIoType::UCHAR: DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback); break; case DcmIoType::CHAR: DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback); break; case DcmIoType::USHORT: DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback); break; case DcmIoType::SHORT: DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback); break; case DcmIoType::UINT: DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback); break; case DcmIoType::INT: DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback); break; case DcmIoType::ULONG: DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback); break; case DcmIoType::LONG: DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback); break; case DcmIoType::FLOAT: DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback); break; case DcmIoType::DOUBLE: DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, correctTilt, callback); break; default: MITK_ERROR << "Found unsupported DICOM pixel type: (enum value) " << io->GetComponentType(); } 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) // transversal 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_NumberOfSlicesApart(1) { } DicomSeriesReader::GantryTiltInformation::GantryTiltInformation( const Point3D& origin1, const Point3D& origin2, const Vector3D& right, const Vector3D& up, unsigned int numberOfSlicesApart) :m_NumberOfSlicesApart(numberOfSlicesApart) { assert(numberOfSlicesApart); // 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 ) */ Vector3D normal = itk::CrossProduct(right, up); Point3D 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; Point3D projectionRight = projectPointOnLine( origin1, origin2, right ); Point3D projectionUp = projectPointOnLine( origin1, origin2, up ); Point3D projectionNormal = projectPointOnLine( origin1, origin2, normal ); m_ShiftRight = (projectionRight - origin2).GetNorm(); m_ShiftUp = (projectionUp - origin2).GetNorm(); m_ShiftNormal = (projectionNormal - origin2).GetNorm(); MITK_DEBUG << " shift normal: " << m_ShiftNormal; MITK_DEBUG << " shift up: " << m_ShiftUp; MITK_DEBUG << " shift right: " << m_ShiftRight; MITK_DEBUG << " tilt angle (rad): " << tanh( m_ShiftUp / m_ShiftNormal ); MITK_DEBUG << " tilt angle (deg): " << tanh( m_ShiftUp / m_ShiftNormal ) * 180.0 / 3.1415926535; } } Point3D DicomSeriesReader::GantryTiltInformation::projectPointOnLine( Point3D p, Point3D lineOrigin, Vector3D lineDirection ) { /** See illustration at http://mo.mathematik.uni-stuttgart.de/inhalt/aussage/aussage472/ vector(lineOrigin,p) = normal * ( innerproduct((p - lineOrigin),normal) / squared-length(normal) ) */ Vector3D lineOriginToP = p - lineOrigin; ScalarType innerProduct = lineOriginToP * lineDirection; ScalarType factor = innerProduct / lineDirection.GetSquaredNorm(); Point3D projection = lineOrigin + factor * lineDirection; return projection; } ScalarType DicomSeriesReader::GantryTiltInformation::GetTiltCorrectedAdditionalSize() const { - return int(m_ShiftUp + 1.0); // to next bigger int - } + // this seems to be a bit too much sometimes, but better too much than cutting off parts of the image + return int(m_ShiftUp + 1.0); // to next bigger int: plus 1, then cut off after point +} ScalarType -DicomSeriesReader::GantryTiltInformation::GetMatrixCoefficientForCorrection() const +DicomSeriesReader::GantryTiltInformation::GetMatrixCoefficientForCorrectionInWorldCoordinates() const { - return - m_ShiftUp / m_ShiftNormal; + // so many mm shifted per slice! + return m_ShiftUp / static_cast(m_NumberOfSlicesApart); } ScalarType DicomSeriesReader::GantryTiltInformation::GetRealZSpacing() const { return m_ShiftNormal / static_cast(m_NumberOfSlicesApart); } bool DicomSeriesReader::GantryTiltInformation::IsSheared() const { return ( m_ShiftRight > 0.001 || m_ShiftUp > 0.001); } bool DicomSeriesReader::GantryTiltInformation::IsRegularGantryTilt() const { return ( m_ShiftRight < 0.001 && m_ShiftUp > 0.001); } std::string DicomSeriesReader::ConstCharStarToString(const char* s) { return s ? std::string(s) : std::string(); } 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, '\\' ) ) { if (dim > 4) break; // otherwise we access invalid array index p[dim++] = atof(coordinate.c_str()); } if (dim != 3) { successful = false; MITK_ERROR << "Reader implementation made wrong assumption on tag (0020,0032). Found " << dim << " instead of 3 values."; } 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, '\\' ) ) { if (dim > 6) break; // otherwise we access invalid array index if (dim<3) { right[dim++] = atof(coordinate.c_str()); } else { up[dim++ - 3] = atof(coordinate.c_str()); } } if (dim != 6) { successful = false; MITK_ERROR << "Reader implementation made wrong assumption on tag (0020,0037). Found " << dim << " instead of 6 values."; } } 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 Vector3D fromFirstToSecondOrigin; fromFirstToSecondOrigin.Fill(0.0); bool fromFirstToSecondOriginInitialized(false); Point3D thisOrigin; 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] ); 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; // Now make sure this direction is along the normal vector of the first slice // If this is NOT the case, then we have a data set with a TILTED GANTRY geometry, // which cannot be loaded into a single mitk::Image at the moment 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() ) { - result.FlagGantryTilt( fromFirstToSecondOrigin ); + result.FlagGantryTilt(); result.AddFileToSortedBlock(*fileIter); // this file is good for current block fileFitsIntoPattern = true; } else { result.AddFileToUnsortedBlock( *fileIter ); // sort away for further analysis fileFitsIntoPattern = false; } } else { 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::UidFileNamesMap DicomSeriesReader::GetSeries(const StringContainer& files, bool groupImagesWithGantryTilt, const StringContainer &restrictions) { return GetSeries(files, true, groupImagesWithGantryTilt, restrictions); } DicomSeriesReader::UidFileNamesMap 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 */ UidFileNamesMap groupsOfSimilarImages; // preliminary result, refined into the final result mapOf3DPlusTBlocks // 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 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 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 ); // 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 // let GDCM scan files if ( !scanner.Scan( files ) ) { MITK_ERROR << "gdcm::Scanner failed when scanning " << files.size() << " input files."; return groupsOfSimilarImages; } // assign files IDs that will separate them for loading into image blocks for (gdcm::Scanner::ConstIterator fileIter = scanner.Begin(); fileIter != scanner.End(); ++fileIter) { //MITK_DEBUG << "Scan file " << fileIter->first << std::endl; if ( std::string(fileIter->first).empty() ) continue; // TODO understand why Scanner has empty string entries // 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) ); groupsOfSimilarImages [ moreUniqueSeriesId ].push_back( fileIter->first ); } // PART II: sort slices spatially for ( UidFileNamesMap::const_iterator groupIter = groupsOfSimilarImages.begin(); groupIter != groupsOfSimilarImages.end(); ++groupIter ) { try { groupsOfSimilarImages[ groupIter->first ] = SortSeriesSlices( groupIter->second ); // sort each slice group spatially } catch(...) { MITK_ERROR << "Catched 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 UidFileNamesMap mapOf3DPlusTBlocks; // final result of this function for ( UidFileNamesMap::const_iterator groupIter = groupsOfSimilarImages.begin(); groupIter != groupsOfSimilarImages.end(); ++groupIter ) { UidFileNamesMap mapOf3DBlocks; // intermediate result for only this group(!) std::map mapOf3DBlockAnalysisResults; StringContainer filesStillToAnalyze = groupIter->second; 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; mapOf3DBlocks[ newGroupUID.str() ] = analysisResult.GetBlockFilenames(); MITK_DEBUG << "Result: sorted 3D group " << newGroupUID.str() << " with " << mapOf3DBlocks[ newGroupUID.str() ].size() << " files"; ++subgroup; filesStillToAnalyze = analysisResult.GetUnsortedFilenames(); // remember what needs further analysis } // end of grouping, now post-process groups // PART IV: attempt to group blocks to 3D+t blocks if requested // inspect entries of mapOf3DBlocks // - 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 // TODO avoid collisions (or prove impossibility) mapOf3DPlusTBlocks.insert( mapOf3DBlocks.begin(), mapOf3DBlocks.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 ( UidFileNamesMap::const_iterator block3DIter = mapOf3DBlocks.begin(); block3DIter != mapOf3DBlocks.end(); ++block3DIter ) { unsigned int numberOfFilesInThisBlock = block3DIter->second.size(); std::string thisBlockKey = block3DIter->first; if (numberOfFilesInPreviousBlock == 0) { numberOfFilesInPreviousBlock = numberOfFilesInThisBlock; mapOf3DPlusTBlocks[thisBlockKey].insert( mapOf3DPlusTBlocks[thisBlockKey].end(), block3DIter->second.begin(), block3DIter->second.end() ); MITK_DEBUG << " 3D+t group " << thisBlockKey << " started"; 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( mapOf3DBlocks[thisBlockKey].front().c_str(), tagImagePositionPatient ), *previous_origin_value = scanner.GetValue( mapOf3DBlocks[previousBlockKey].front().c_str(), tagImagePositionPatient ), *destination_value = scanner.GetValue( mapOf3DBlocks[thisBlockKey].back().c_str(), tagImagePositionPatient ), *previous_destination_value = scanner.GetValue( mapOf3DBlocks[previousBlockKey].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 mapOf3DPlusTBlocks[previousBlockKey].insert( mapOf3DPlusTBlocks[previousBlockKey].end(), block3DIter->second.begin(), block3DIter->second.end() ); MITK_DEBUG << " --> group enhanced with another timestep"; } else { // start a new block mapOf3DPlusTBlocks[thisBlockKey].insert( mapOf3DPlusTBlocks[thisBlockKey].end(), block3DIter->second.begin(), block3DIter->second.end() ); MITK_DEBUG << " ==> group closed with " << mapOf3DPlusTBlocks[previousBlockKey].size() / numberOfFilesInPreviousBlock << " time steps"; previousBlockKey = thisBlockKey; MITK_DEBUG << " 3D+t group " << thisBlockKey << " started"; } } numberOfFilesInPreviousBlock = numberOfFilesInThisBlock; } } } MITK_DEBUG << "================================================================================"; MITK_DEBUG << "Summary: "; for ( UidFileNamesMap::const_iterator groupIter = mapOf3DPlusTBlocks.begin(); groupIter != mapOf3DPlusTBlocks.end(); ++groupIter ) { MITK_DEBUG << " Image volume " << groupIter->first << " with " << groupIter->second.size() << " files"; } MITK_DEBUG << "Done. "; MITK_DEBUG << "================================================================================"; return mapOf3DPlusTBlocks; } DicomSeriesReader::UidFileNamesMap 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& e) { 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 tagSliceThickness(0x0018, 0x0050); // slice thickness const gdcm::Tag tagNumberOfRows(0x0028, 0x0010); // number rows const gdcm::Tag tagNumberOfColumns(0x0028, 0x0011); // number cols std::string constructedID; try { constructedID = tagValueMap[ tagSeriesInstanceUID ]; } catch (std::exception& e) { MITK_ERROR << "CreateMoreUniqueSeriesIdentifier() could not access series instance UID. Something is seriously wrong with this image."; MITK_ERROR << "Error from exception: " << e.what(); } constructedID += CreateSeriesIdentifierPart( tagValueMap, tagNumberOfRows ); constructedID += CreateSeriesIdentifierPart( tagValueMap, tagNumberOfColumns ); constructedID += CreateSeriesIdentifierPart( tagValueMap, tagPixelSpacing ); constructedID += CreateSeriesIdentifierPart( tagValueMap, tagSliceThickness ); constructedID += CreateSeriesIdentifierPart( tagValueMap, tagImageOrientation ); 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) { UidFileNamesMap allSeries = GetSeries(dir, groupImagesWithGantryTilt, restrictions); StringContainer resultingFileList; for ( UidFileNamesMap::const_iterator idIter = allSeries.begin(); idIter != allSeries.end(); ++idIter ) { if ( idIter->first.find( series_uid ) == 0 ) // this ID starts with given series_uid { resultingFileList.insert( resultingFileList.end(), idIter->second.begin(), idIter->second.end() ); // append } } return resultingFileList; } DicomSeriesReader::StringContainer DicomSeriesReader::SortSeriesSlices(const StringContainer &unsortedFilenames) { gdcm::Sorter sorter; sorter.SetSortFunction(DicomSeriesReader::GdcmSortFunction); try { sorter.Sort(unsortedFilenames); return sorter.GetFilenames(); } catch(std::logic_error&) { MITK_WARN << "Sorting error. Leaving series unsorted."; return unsortedFilenames; } } bool DicomSeriesReader::GdcmSortFunction(const gdcm::DataSet &ds1, const gdcm::DataSet &ds2) { // make sure we habe Image Position and Orientation if ( ! ( ds1.FindDataElement(gdcm::Tag(0x0020,0x0032)) && ds1.FindDataElement(gdcm::Tag(0x0020,0x0037)) && ds2.FindDataElement(gdcm::Tag(0x0020,0x0032)) && ds2.FindDataElement(gdcm::Tag(0x0020,0x0037)) ) ) { MITK_WARN << "Dicom images are missing attributes for a meaningful sorting."; throw std::logic_error("Dicom images are missing attributes for a meaningful sorting."); } 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); if (image_orientation1 != image_orientation2) { 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; for (unsigned char i = 0u; i < 3u; ++i) { dist1 += normal[i] * image_pos1[i]; dist2 += normal[i] * image_pos2[i]; } if ( fabs(dist1 - dist2) < mitk::eps) { gdcm::Attribute<0x0008,0x0032> acq_time1; // Acquisition time (may be missing, so we check existence first) gdcm::Attribute<0x0008,0x0032> acq_time2; gdcm::Attribute<0x0020,0x0012> acq_number1; // Acquisition number (may also be missing, so we check existence first) gdcm::Attribute<0x0020,0x0012> acq_number2; if (ds1.FindDataElement(gdcm::Tag(0x0008,0x0032)) && ds2.FindDataElement(gdcm::Tag(0x0008,0x0032))) { acq_time1.Set(ds1); acq_time2.Set(ds2); return acq_time1 < acq_time2; } else if (ds1.FindDataElement(gdcm::Tag(0x0020,0x0012)) && ds2.FindDataElement(gdcm::Tag(0x0020,0x0012))) { acq_number1.Set(ds1); acq_number2.Set(ds2); return acq_number1 < acq_number2; } else { return true; } } else { // default: compare position return dist1 < dist2; } } 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, Image *image) { std::list imageBlock; imageBlock.push_back(filenames); CopyMetaDataToImageProperties(imageBlock, tagValueMappings_, io, image); } void DicomSeriesReader::CopyMetaDataToImageProperties( std::list imageBlock, const gdcm::Scanner::MappingType& tagValueMappings_, DcmIoType* io, 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) { propertyKeySliceLocation.append(".t" + timeStep); propertyKeyInstanceNumber.append(".t" + timeStep); propertyKeySOPInstanceNumber.append(".t" + timeStep); } 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; } } } // end namespace mitk #include diff --git a/Core/Code/IO/mitkDicomSeriesReader.h b/Core/Code/IO/mitkDicomSeriesReader.h index f96f699cfc..d7484fb62d 100644 --- a/Core/Code/IO/mitkDicomSeriesReader.h +++ b/Core/Code/IO/mitkDicomSeriesReader.h @@ -1,612 +1,687 @@ /*=================================================================== 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 #include #include #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_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 correcness of generated slice positions! \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 for SOP Classes CT Image Storage and MR Image Storage (NOT for the "Enhanced" variants containing multi-frame images) - 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 \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 - Images from tilted CT gantries CAN ONLY be loaded as a series of single-slice images, since mitk::Image or the accompanying mapper are not (yet?) capable of representing such geometries - Secondary Capture images are expected to have the (0018,2010) tag describing the pixel spacing. If only the (0028,0030) tag is set, the spacing will be misinterpreted as (1,1) \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 identifiers similar to DICOM UIDs, each associated to a sorted list of file names. 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 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::UidFileNamesMap 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,0037) Image Orientation - (0028,0030) Pixel Spacing - (0018,0050) Slice Thickness - (0028,0010) Number Of Rows - (0028,0011) Number Of Columns - (0020,000e) Series Instance UID : could be argued about, might be dropped in the future (optionally) \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, (0008,0032) Acquisition Time is used as a backup criterion (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. 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 slices of a series, both described as a pair (origin/orientation): - we calculate if the second origin is on a line along the normal of the first 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. \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 For grouped lists of filenames, assigned an ID each. */ typedef std::map UidFileNamesMap; /** \brief Interface for the progress callback. */ typedef void (*UpdateCallBackMethod)(float); /** \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 UidFileNamesMap 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 UidFileNamesMap 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 UidFileNamesMap 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); /** \brief See LoadDicomSeries! Just a slightly different interface. */ static bool LoadDicomSeries(const StringContainer &filenames, DataNode &node, bool sort = true, bool load4D = true, bool correctGantryTilt = true, UpdateCallBackMethod callback = 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(); - StringContainer GetUnsortedFilenames(); + /** + \brief Remaining files, which could not be grouped. + */ + StringContainer GetUnsortedFilenames(); + + /** + \brief Wheter or not the grouped result contain a gantry tilt. + */ bool ContainsGantryTilt(); - Vector3D GetInterSliceOffset(); + /** + \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 FlagGantryTilt(Vector3D interSliceOffset); - + + /** + \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; - - Vector3D m_InterSliceOffset; }; + /** + \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. + + All calculations are done in the constructor, results can then + be read via the remaining methods. + */ class GantryTiltInformation { public: + /** + \brief Just so we can create empty instances for assigning results later. + */ GantryTiltInformation(); - + + /** + \brief THE constructor, which does all the calculations. + + See code comments for explanation. + */ GantryTiltInformation( const Point3D& origin1, const Point3D& origin2, const Vector3D& right, const Vector3D& up, unsigned int numberOfSlicesApart); - + + /** + \brief Whether the slices were sheared. + */ bool IsSheared() const; + + /** + \brief Whether the shearing is a gantry tilt or more complicated. + */ bool IsRegularGantryTilt() const; - ScalarType GetMatrixCoefficientForCorrection() const; + /** + \brief The offset distance in Y direction for each slice (describes the tilt result). + */ + ScalarType GetMatrixCoefficientForCorrectionInWorldCoordinates() const; + + + /** + \brief The z / inter-slice spacing. Needed to correct ImageSeriesReader's result. + */ ScalarType GetRealZSpacing() const; + /** + \brief The shift between first and last slice in mm. + + Needed to resize an orthogonal image volume. + */ ScalarType GetTiltCorrectedAdditionalSize() const; protected: + /** + \brief Projection of point p onto line through lineOrigin in direction of lineDirection. + */ Point3D projectPointOnLine( Point3D p, Point3D lineOrigin, Vector3D lineDirection ); ScalarType m_ShiftUp; ScalarType m_ShiftRight; ScalarType m_ShiftNormal; 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_); static std::string ConstCharStarToString(const char* s); static Point3D DICOMStringToPoint3D(const std::string& s, bool& successful); static void DICOMStringToOrientationVectors(const std::string& s, Vector3D& right, Vector3D& up, bool& successful); template static typename ImageType::Pointer 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); 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; }; /** \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); /** \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, CallbackCommand* command = NULL); /** \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, Image* image); static void CopyMetaDataToImageProperties( std::list imageBlock, const gdcm::Scanner::MappingType& tagValueMappings_, DcmIoType* io, 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 084c23bed8..80c95e54c9 100644 --- a/Core/Code/IO/mitkDicomSeriesReader.txx +++ b/Core/Code/IO/mitkDicomSeriesReader.txx @@ -1,455 +1,489 @@ /*=================================================================== 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 namespace mitk { template void DicomSeriesReader::LoadDicom(const StringContainer &filenames, DataNode &node, bool sort, bool load4D, bool correctTilt, UpdateCallBackMethod callback) { 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); const gdcm::Tag tagImagePositionPatient(0x0020,0x0032); // Image Position (Patient) const gdcm::Tag tagImageOrientation(0x0020, 0x0037); // Image Orientation try { mitk::Image::Pointer image = mitk::Image::New(); 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()) ) { 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(); 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(); MITK_DEBUG << "** Loading now: shear? " << tiltInfo.IsSheared(); MITK_DEBUG << "** Loading now: normal tilt? " << tiltInfo.IsRegularGantryTilt(); MITK_DEBUG << "** Loading now: perform tilt correction? " << correctTilt; } else { correctTilt = false; // we CANNOT do that } if (volume_count == 1 || !canLoadAs4D || !load4D) { image = LoadDICOMByITK( imageBlocks.front(), correctTilt, tiltInfo, command ); // load first 3D block initialize_node = true; } else if (volume_count > 1) { // 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; 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); gdcm::Scanner scanner; ScanForSliceInformation(filenames, scanner); DicomSeriesReader::CopyMetaDataToImageProperties( imageBlocks, scanner.GetMappings(), io, image); MITK_DEBUG << "Volume dimension: [" << image->GetDimension(0) << ", " << image->GetDimension(1) << ", " << image->GetDimension(2) << ", " << image->GetDimension(3) << "]"; #if (GDCM_MAJOR_VERSION == 2) && (GDCM_MINOR_VERSION < 1) && (GDCM_BUILD_VERSION < 15) // workaround for a GDCM 2 bug until version 2.0.15: // GDCM read spacing vector wrongly. Instead of "row spacing, column spacing", it misinterprets the DICOM tag as "column spacing, row spacing". // this is undone here, until we use a GDCM that has this issue fixed. // From the commit comments, GDCM 2.0.15 fixed the spacing interpretation with bug 2901181 // http://sourceforge.net/tracker/index.php?func=detail&aid=2901181&group_id=137895&atid=739587 Vector3D correctedImageSpacing = image->GetGeometry()->GetSpacing(); std::swap( correctedImageSpacing[0], correctedImageSpacing[1] ); image->GetGeometry()->SetSpacing( correctedImageSpacing ); #endif MITK_DEBUG << "Volume spacing: [" << image->GetGeometry()->GetSpacing()[0] << ", " << image->GetGeometry()->GetSpacing()[1] << ", " << image->GetGeometry()->GetSpacing()[2] << "]"; for (std::list::iterator df_it = ++imageBlocks.begin(); df_it != imageBlocks.end(); ++df_it) { reader->SetFileNames(*df_it); reader->Update(); readVolume = reader->GetOutput(); if (correctTilt) { readVolume = InPlaceFixUpTiltedGeometry( reader->GetOutput(), tiltInfo ); } image->SetImportVolume(readVolume->GetBufferPointer(), act_volume++); } 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); } } catch (std::exception& e) { // reset locale then throw up setlocale(LC_NUMERIC, previousCLocale); std::cin.imbue(previousCppLocale); throw e; } } template Image::Pointer DicomSeriesReader::LoadDICOMByITK( const StringContainer& filenames, bool correctTilt, const GantryTiltInformation& tiltInfo, CallbackCommand* command ) { /******** Normal Case, 3D (also for GDCM < 2 usable) ***************/ mitk::Image::Pointer image = mitk::Image::New(); 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); } 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()); gdcm::Scanner scanner; ScanForSliceInformation(filenames, scanner); DicomSeriesReader::CopyMetaDataToImageProperties( filenames, scanner.GetMappings(), io, image); MITK_DEBUG << "Volume dimension: [" << image->GetDimension(0) << ", " << image->GetDimension(1) << ", " << image->GetDimension(2) << "]"; #if (GDCM_MAJOR_VERSION == 2) && (GDCM_MINOR_VERSION < 1) && (GDCM_BUILD_VERSION < 15) // workaround for a GDCM 2 bug until version 2.0.15: // GDCM read spacing vector wrongly. Instead of "row spacing, column spacing", it misinterprets the DICOM tag as "column spacing, row spacing". // this is undone here, until we use a GDCM that has this issue fixed. // From the commit comments, GDCM 2.0.15 fixed the spacing interpretation with bug 2901181 // http://sourceforge.net/tracker/index.php?func=detail&aid=2901181&group_id=137895&atid=739587 Vector3D correctedImageSpacing = image->GetGeometry()->GetSpacing(); std::swap( correctedImageSpacing[0], correctedImageSpacing[1] ); image->GetGeometry()->SetSpacing( correctedImageSpacing ); #endif 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 ippTag(0x0020,0x0032); //Image position (Patient) scanner.AddTag(ippTag); const gdcm::Tag tagImageOrientation(0x0020,0x0037); //Image orientation scanner.AddTag(tagImageOrientation); // TODO what if tags don't exist? 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 ); scanner.Scan(filenames); // make available image position 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 const gdcm::Tag ippTag(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(ippTag) == tagToValueMap.end()) { continue; } std::string position = tagToValueMap.find(ippTag)->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 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 + Anybody who does this in a simpler way: don't forget to write up how and why your solution works */ typedef itk::AffineTransform< double, ImageType::ImageDimension > TransformType; typename TransformType::Pointer transformShear = TransformType::New(); + /** + What I think should be done here: + + - 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 should bring 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 probably would need the shift factor in index coordinates but we use mm, which appears strange, see below + - we transform the image volume back to its actual position + - presumably back 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 + + Here comes the unsolved PROBLEM: + - WHY is it correct to apply the shift factor in millimeters, when we assume we are in a index coordinate system? + - or why is it not millimeters but index coordinates? + - or what else is the wrong assumption here? + - This is a serious question to anybody very firm in transforms, ITK, etc. + - this is also a chocolate reward question. Provide a good explanation with corrected code as a git commit and get it. + */ + + // ScalarType factor = tiltInfo.GetMatrixCoefficientForCorrectionInWorldCoordinates() / input->GetSpacing()[1]; + ScalarType factor = tiltInfo.GetMatrixCoefficientForCorrectionInWorldCoordinates(); // row 1, column 2 corrects shear in parallel to Y axis, proportional to distance in Z direction - transformShear->Shear( 1, 2, tiltInfo.GetMatrixCoefficientForCorrection() ); + transformShear->Shear( 1, 2, factor ); typename TransformType::Pointer imageIndexToWorld = TransformType::New(); imageIndexToWorld->SetOffset( input->GetOrigin().GetVectorFromOrigin() ); imageIndexToWorld->SetMatrix( input->GetDirection() ); 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 ); - resampler->SetDefaultPixelValue( std::numeric_limits::min() ); // TODO minimum meaningful modality value? + /* + 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. + */ + resampler->SetDefaultPixelValue( std::numeric_limits::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 typename ImageType::SizeType largerSize = resampler->GetSize(); // now the resampler already holds the input image's size. largerSize[1] += tiltInfo.GetTiltCorrectedAdditionalSize(); resampler->SetSize( largerSize ); 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