diff --git a/Core/Code/IO/mitkDicomSeriesReader.cpp b/Core/Code/IO/mitkDicomSeriesReader.cpp index d3c76c0545..909224b4a3 100644 --- a/Core/Code/IO/mitkDicomSeriesReader.cpp +++ b/Core/Code/IO/mitkDicomSeriesReader.cpp @@ -1,760 +1,903 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date$ Version: $Revision$ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ //#define MBILOG_ENABLE_DEBUG #include #include #if GDCM_MAJOR_VERSION >= 2 #include #include #include #include #include #endif namespace mitk { typedef itk::GDCMSeriesFileNames DcmFileNamesGeneratorType; DataNode::Pointer DicomSeriesReader::LoadDicomSeries(const StringContainer &filenames, bool sort, bool check_4d, UpdateCallBackMethod callback) { DataNode::Pointer node = DataNode::New(); if (DicomSeriesReader::LoadDicomSeries(filenames, *node, sort, check_4d, callback)) { if( filenames.empty() ) { return NULL; } return node; } else { return NULL; } } bool DicomSeriesReader::LoadDicomSeries(const StringContainer &filenames, DataNode &node, bool sort, bool check_4d, 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, callback); return true; case DcmIoType::CHAR: DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, callback); return true; case DcmIoType::USHORT: DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, callback); return true; case DcmIoType::SHORT: DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, callback); return true; case DcmIoType::UINT: DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, callback); return true; case DcmIoType::INT: DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, callback); return true; case DcmIoType::ULONG: DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, callback); return true; case DcmIoType::LONG: DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, callback); return true; case DcmIoType::FLOAT: DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, callback); return true; case DcmIoType::DOUBLE: DicomSeriesReader::LoadDicom(filenames, node, sort, check_4d, callback); return true; default: MITK_ERROR << "Found unsupported DICOM pixel type: (enum value) " << io->GetComponentType(); } } } 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()); } #if GDCM_MAJOR_VERSION >= 2 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 } #endif DicomSeriesReader::TwoStringContainers DicomSeriesReader::AnalyzeFileForITKImageSeriesReaderSpacingAssumption( const StringContainer& files, const gdcm::Scanner::MappingType& tagValueMappings_) { // result.first = files that fit ITK's assumption // result.second = files that do not fit TwoStringContainers 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); - Point3D thisOrigin; thisOrigin.Fill(0.0); - Point3D lastOrigin; lastOrigin.Fill(0.0); + bool fromFirstToSecondOriginInitialized(false); + Point3D thisOrigin; + Point3D lastOrigin; + Point3D lastDifferentOrigin; + bool lastOriginInitialized(false); MITK_DEBUG << "Analyzing files for z-spacing assumption of ITK's ImageSeriesReader "; unsigned int fileIndex(0); for (StringContainer::const_iterator fileIter = files.begin(); fileIter != files.end(); ++fileIter, ++fileIndex) { // Read tag value into point3D. PLEASE replace this by appropriate GDCM code if you figure out how to do that std::string thisOriginString = tagValueMappings[fileIter->c_str()][tagImagePositionPatient]; std::istringstream originReader(thisOriginString); std::string coordinate; unsigned int dim(0); while( std::getline( originReader, coordinate, '\\' ) ) { thisOrigin[dim++] = atof(coordinate.c_str()); } if (dim != 3) { MITK_ERROR << "Reader implementation made wrong assumption on tag (0020,0032). Found " << dim << "instead of 3 values."; } MITK_DEBUG << " " << fileIndex << " " << *fileIter << " at " << thisOriginString << "(" << thisOrigin[0] << "," << thisOrigin[1] << "," << thisOrigin[2] << ")"; - // fill in vector value as soon as possible - if (fileIndex == 1) + if ( lastOriginInitialized && (thisOrigin == lastOrigin) ) { - fromFirstToSecondOrigin = thisOrigin - lastOrigin; + MITK_DEBUG << "Sort away " << *fileIter << " for separate time step"; + // TODO sort away this slice, we already have one occupying this position + result.second.push_back( *fileIter ); } - else if (fileIndex > 1) + else { - Point3D assumedOrigin = lastOrigin + fromFirstToSecondOrigin; + // fill in vector value as soon as possible + if (!fromFirstToSecondOriginInitialized && lastOriginInitialized) // calculate vector as soon as possible when we get a new position + { + fromFirstToSecondOrigin = thisOrigin - lastDifferentOrigin; + fromFirstToSecondOriginInitialized = true; - Vector3D originError = assumedOrigin - thisOrigin; - double norm = originError.GetNorm(); + // 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 - if (norm > 3 * mitk::sqrteps) - { - MITK_WARN << "File " << *fileIter << " breaks the inter-slice distance pattern (diff = " - << norm << ", allowed " - << 3 * (mitk::sqrteps)<< ")."; - MITK_WARN << "Expected position (" << assumedOrigin[0] << "," - << assumedOrigin[1] << "," - << assumedOrigin[2] << "), got position (" - << thisOrigin[0] << "," - << thisOrigin[1] << "," - << thisOrigin[2] << ")"; - - // 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 + // Again ugly code to read tag Image Orientation into two Vectors + Vector3D right; right.Fill(0.0); + Vector3D up; right.Fill(0.0); // might be down as, but it is just a name at this point + std::string thisOrientationString = tagValueMappings[fileIter->c_str()][tagImageOrientation]; - result.first.assign( files.begin(), fileIter ); - result.second.assign( fileIter, files.end() ); + std::istringstream orientationReader(thisOrientationString); + std::string coordinate; + unsigned int dim(0); + while( std::getline( orientationReader, coordinate, '\\' ) ) + { + if (dim<3) + { + right[dim++] = atof(coordinate.c_str()); + } + else + { + up[dim++ - 3] = atof(coordinate.c_str()); + } + } + + if (dim != 6) + { + MITK_ERROR << "Reader implementation made wrong assumption on tag (0020,0037). Found " << dim << "instead of 6 values."; + } + + MITK_DEBUG << "Tilt check: right vector (" << right[0] << "," << right[1] << "," << right[2] << "), " + "up vector (" << up[0] << "," << up[1] << "," << up[2] << ")"; + + /* + Determine if line (thisOrigin + l * normal) contains lastDifferentOrigin. + + Done by calculating the distance of lastDifferentOrigin from line (thisOrigin + l *normal) + + E.g. http://mathworld.wolfram.com/Point-LineDistance3-Dimensional.html + + squared distance = | (pointAlongNormal - thisOrign) x (thisOrigin - lastDifferentOrigin) | ^ 2 + / + |pointAlongNormal - thisOrigin| ^ 2 + + ( x meaning the cross product ) + */ + + Vector3D normal = itk::CrossProduct(right, up); + Point3D pointAlongNormal = thisOrigin + normal; + + double numerator = itk::CrossProduct( pointAlongNormal - thisOrigin , thisOrigin - lastDifferentOrigin ).GetSquaredNorm(); + double denominator = (pointAlongNormal - thisOrigin).GetSquaredNorm(); + + double distance = sqrt(numerator / denominator); - return result; // stop processing with first split + if (distance > 0.001) + { + MITK_WARN << "Series seems to contain a tilted geometry. Will load as many single slices."; + MITK_WARN << "Distance of expected origin from actual origin: " << distance; + + result.first.assign( files.begin(), fileIter ); + result.second.insert( result.second.end(), fileIter, files.end() ); + + return result; // stop processing with first split + } + } + else if (fromFirstToSecondOriginInitialized) // we already know the offset between slices + { + Point3D assumedOrigin = lastDifferentOrigin + fromFirstToSecondOrigin; + + Vector3D originError = assumedOrigin - thisOrigin; + double norm = originError.GetNorm(); + + if (norm > 3 * mitk::sqrteps) + { + MITK_WARN << "File " << *fileIter << " breaks the inter-slice distance pattern (diff = " + << norm << ", allowed " + << 3 * (mitk::sqrteps)<< ")."; + MITK_WARN << "Expected position (" << assumedOrigin[0] << "," + << assumedOrigin[1] << "," + << assumedOrigin[2] << "), got position (" + << thisOrigin[0] << "," + << thisOrigin[1] << "," + << thisOrigin[2] << ")"; + + // 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 + + result.first.assign( files.begin(), fileIter ); + result.second.insert( result.second.end(), fileIter, files.end() ); + + return result; // stop processing with first split + } } + + result.first.push_back(*fileIter); } - lastOrigin = thisOrigin; + // recored current origin for reference in later iterations + if ( !lastOriginInitialized || thisOrigin != lastOrigin ) + { + lastDifferentOrigin = thisOrigin; + } + + lastOrigin = thisOrigin; + lastOriginInitialized = true; } - // default: all files match our expectations: return unchanged list - result.first = files; - return result; } - DicomSeriesReader::UidFileNamesMap DicomSeriesReader::GetSeries(const StringContainer& files, const StringContainer &restrictions) +{ + return GetSeries(files, true, restrictions); +} + +DicomSeriesReader::UidFileNamesMap +DicomSeriesReader::GetSeries(const StringContainer& files, bool sortTo3DPlust, const StringContainer &restrictions) { UidFileNamesMap map; // result variable /** 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 */ #if GDCM_MAJOR_VERSION < 2 // old GDCM: let itk::GDCMSeriesFileNames do the sorting DcmFileNamesGeneratorType::Pointer name_generator = DcmFileNamesGeneratorType::New(); name_generator->SetUseSeriesDetails(true); name_generator->AddSeriesRestriction("0020|0037"); // image orientation (patient) name_generator->AddSeriesRestriction("0028|0030"); // pixel spacing (x,y) name_generator->SetLoadSequences(false); // could speed up reading, and we don't use sequences anyway name_generator->SetLoadPrivateTags(false); for(StringContainer::const_iterator it = restrictions.begin(); it != restrictions.end(); ++it) { name_generator->AddSeriesRestriction(*it); } name_generator->SetDirectory(dir.c_str()); const StringContainer& series_uids = name_generator->GetSeriesUIDs(); for(StringContainer::const_iterator it = series_uids.begin(); it != series_uids.end(); ++it) { const std::string& uid = *it; map[uid] = name_generator->GetFileNames(uid); } #else // new GDCM: 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 map; } // 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 << "Read 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) ); map [ moreUniqueSeriesId ].push_back( fileIter->first ); } // PART II: analyze pre-sorted images for valid blocks, // separate into multiple blocks if necessary. // // Analysis performs the following steps: // * sort slices spatially // * 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 for ( UidFileNamesMap::const_iterator groupIter = map.begin(); groupIter != map.end(); ++groupIter ) { // sort each slice group spatially map[ groupIter->first ] = SortSeriesSlices( groupIter->second ); } - UidFileNamesMap enhancedMap; // additional results generated by this analysis step + UidFileNamesMap mapOf3DPlusTBlocks; // sorting results generated by this analysis step for ( UidFileNamesMap::const_iterator groupIter = map.begin(); groupIter != map.end(); ++groupIter ) { + UidFileNamesMap mapOf3DBlocks; // intermediate result for only this group(!) 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 { TwoStringContainers analysisResult = AnalyzeFileForITKImageSeriesReaderSpacingAssumption( filesStillToAnalyze, scanner.GetMappings() ); - if (subgroup > 0) // record what has already been sorted out - { - // enhance the UID for additional groups - std::stringstream newGroupUID; - newGroupUID << groupUID << '.' << subgroup; - enhancedMap[ newGroupUID.str() ] = analysisResult.first; - MITK_INFO << "Sorted group " << newGroupUID.str() << " with " << enhancedMap[ newGroupUID.str() ].size() << " files"; - } - else - { - map[ groupUID ] = analysisResult.first; - MITK_INFO << "Sorted group " << groupUID << " with " << map[ groupUID ].size() << " files"; - } + // enhance the UID for additional groups + std::stringstream newGroupUID; + newGroupUID << groupUID << '.' << subgroup; + mapOf3DBlocks[ newGroupUID.str() ] = analysisResult.first; + MITK_INFO << "Sorted 3D group " << newGroupUID.str() << " with " << mapOf3DBlocks[ newGroupUID.str() ].size() << " files"; + ++subgroup; filesStillToAnalyze = analysisResult.second; // remember what needs further analysis } - } - - // add newly generated results TODO avoid collisions (or prove impossibility) - map.insert( enhancedMap.begin(), enhancedMap.end() ); - /* TODO at this point: + // end of grouping, now post-process groups + + // PART III: 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 - sort slices spatially SortSeriesSlices - sort into blocks of equal spacing + // 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 III") - */ + 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_INFO << "3D+t group " << thisBlockKey << " started"; + previousBlockKey = thisBlockKey; + } + else + { + // check whether this and the previous block share a comon origin + std::string thisOriginString = scanner.GetValue( mapOf3DBlocks[thisBlockKey].front().c_str(), tagImagePositionPatient ); + std::string previousOriginString = scanner.GetValue( mapOf3DBlocks[previousBlockKey].front().c_str(), tagImagePositionPatient ); + // TODO should be safe, but a little try/catch or other error handling wouldn't hurt + + bool identicalOrigins(thisOriginString == previousOriginString); + + if (identicalOrigins && (numberOfFilesInPreviousBlock == numberOfFilesInThisBlock)) + { + // group with previous block + mapOf3DPlusTBlocks[previousBlockKey].insert( mapOf3DPlusTBlocks[previousBlockKey].end(), + block3DIter->second.begin(), + block3DIter->second.end() ); + MITK_INFO << "3D+t group " << previousBlockKey << " enhanced with another timestep"; + } + else + { + // start a new block + mapOf3DPlusTBlocks[thisBlockKey].insert( mapOf3DPlusTBlocks[thisBlockKey].end(), + block3DIter->second.begin(), + block3DIter->second.end() ); + MITK_INFO << "3D+t group " << thisBlockKey << " started"; + previousBlockKey = thisBlockKey; + } + } + + numberOfFilesInPreviousBlock = numberOfFilesInThisBlock; + } + } + } #endif for ( UidFileNamesMap::const_iterator groupIter = map.begin(); groupIter != map.end(); ++groupIter ) { MITK_DEBUG << "Slice group " << groupIter->first << " with " << groupIter->second.size() << " files"; } - return map; + return mapOf3DPlusTBlocks; } DicomSeriesReader::UidFileNamesMap DicomSeriesReader::GetSeries(const std::string &dir, const StringContainer &restrictions) { gdcm::Directory directoryLister; directoryLister.Load( dir.c_str(), false ); // non-recursive return GetSeries(directoryLister.GetFilenames(), restrictions); } #if GDCM_MAJOR_VERSION >= 2 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; } #endif DicomSeriesReader::StringContainer DicomSeriesReader::GetSeries(const std::string &dir, const std::string &series_uid, const StringContainer &restrictions) { UidFileNamesMap allSeries = GetSeries(dir, 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) { #if GDCM_MAJOR_VERSION >= 2 gdcm::Sorter sorter; sorter.SetSortFunction(DicomSeriesReader::GdcmSortFunction); sorter.Sort(unsortedFilenames); return sorter.GetFilenames(); #else return unsortedFilenames; #endif } #if GDCM_MAJOR_VERSION >= 2 bool DicomSeriesReader::GdcmSortFunction(const gdcm::DataSet &ds1, const gdcm::DataSet &ds2) { 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; if (ds1.FindDataElement(gdcm::Tag(0x0008,0x0032))) acq_time1.Set(ds1); if (ds2.FindDataElement(gdcm::Tag(0x0008,0x0032))) acq_time2.Set(ds2); + // TODO this could lead to comparison of unset times (does Attribute initialize to good defaults?) // exception: same position: compare by acquisition time return acq_time1 < acq_time2; } else { // default: compare position return dist1 < dist2; } } #endif std::string DicomSeriesReader::GetConfigurationString() { std::stringstream configuration; configuration << "MITK_USE_GDCMIO: "; #ifdef MITK_USE_GDCMIO configuration << "true"; #else configuration << "false"; #endif configuration << "\n"; configuration << "GDCM_VERSION: "; #ifdef GDCM_MAJOR_VERSION configuration << GDCM_VERSION; #endif //configuration << "\n"; return configuration.str(); } } #include diff --git a/Core/Code/IO/mitkDicomSeriesReader.h b/Core/Code/IO/mitkDicomSeriesReader.h index ab003ebbfa..96cfc4f018 100644 --- a/Core/Code/IO/mitkDicomSeriesReader.h +++ b/Core/Code/IO/mitkDicomSeriesReader.h @@ -1,397 +1,423 @@ /*========================================================================= Program: Medical Imaging & Interaction Toolkit Language: C++ Date: $Date$ Version: $Revision$ Copyright (c) German Cancer Research Center, Division of Medical and Biological Informatics. All rights reserved. See MITKCopyright.txt or http://www.mitk.org/copyright.html for details. This software is distributed WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the above copyright notices for more information. =========================================================================*/ #ifndef mitkDicomSeriesReader_h #define mitkDicomSeriesReader_h #include "mitkDataNode.h" #include "mitkConfig.h" #ifdef MITK_USE_GDCMIO #include #else #include #endif #include #include #include #if GDCM_MAJOR_VERSION >= 2 #include #include #include #include #include #include #endif namespace mitk { /** \brief Central DICOM image loading class for MITK. \section DicomSeriesReader_usage Usage The DICOM series is specified through an enumeration of all files contained in it. A method for the creation of file name enumerations of all series contained in a particular directory is also provided. After loading the series, the generated DataNode is not named. The caller is responsible for given the node a name. Example: \code // only directory known here: /home/me/dicom DicomSeriesReader::UidFileNamesMap allImageBlocks = DicomSeriesReader::GetSeries("/home/me/dicom/"); // file now divided into groups of identical image size, orientation, spacing, etc. // each of these lists should be loadable as an mitk::Image. // optional step: sorting (this is also done by LoadDicomSeries if not explicitly suppressed) DicomSeriesReader::StringContainer seriesToLoad = allImageBlocks[...]; // decide for yourself DicomSeriesReader::StringContainer oneBlockSorted = DicomSeriesReader::SortSeriesSlices( seriesToLoad ); // final step: load into DataNode (can result in 3D+t image) DataNode::Pointer node = DicomSeriesReader::LoadDicomSeries( oneBlockSorted ); // would sort again, but we could set the sort flag to false Image::Pointer image = dynamic_cast( node->GetData() ); \endcode TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO Necessary changes (bug 7285) * GetSeries should by default * sort slices * determine if slice groups contain gaps (this is only done fast if slices are sorted previously) of if slices are not equally spaced * if gap, split into two parts * if non-continous: split into continuous groups * this is a little tricky: create (all?) possible sortings, then determine the one with minimum number of groups * If somebody demands it, this check could be explicitly disabled by a flag * The SortSeriesSlices method can be removed (take care of existing callers), since this is done by GetSeries * LoadDicomSeries changes * remove parameter "bool sort" (now expected to be done previously) * handle titled geometries! * ImageSeriesReader will definitely calculate a wrong spacing * either we load each slice separately * or ...? + + + 3D+t handling: + * the "block sorting" in GetSeries will separate images that occupy the identical position in space + * this process will keep the initial sorting that organized slices by 1. position 2. acquisition time + * i.e. after block sorting, there are multiple blocks that cover identical positions and + have identical properties regarding orientation, dimensions, spacing etc. + * depending on arguments passed to GetSeries, it will attempt to group as many image blocks as possible to + a single, large 3D+t image group (default behavior is to do this). + * LoadDicomSeries expects(!) that the filenames passed either form a valid 3D block or a 3D+t block + as sorted by GetSeries() + + * TODO later, when dicom tags are loaded for each mitk::Image slice, there should be a verification + step for correct positions (optional if time-consuming) + TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO TODO \section DicomSeriesReader_sorting Sorting into 3D+t blocks TODO: describe strategy of sorting images into blocks of 3D image stacks \section DicomSeriesReader_limitations Assumptions and limitations The class is mostly updated for working properly with GDCM 2.0.14 and MITK_USE_GDCMIO ON. The other version are kept here for compatibility and may be removed in the future. \b Assumptions - expected to work for CT Image and MR Image - special treatment for a certain type of Philips 3D ultrasound (recogized by tag 3001,0010 set to "Philips3D") \b Limitations - 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_internal Internal structure only valid for gdcm 2.0.14+ map GetSeries(dir) - scan for DICOM files in dir (gdcm::Scanner) - read Series Instance UID, image orientation, pixel spacing, slice thickness, number rows, number cols - sort all files into groups of files that match in mentioned attributes - resulting groups are named as "SeriesInstanceUID + UID like string made up of attribute values" StringContainer SortSeriesSlices(StringContainer) - use gdcm::Sorter to sort spatially - sort by 1. distance along slice normal direction 2. acquisition time (if present) 3. by input order list SortIntoBlocksFor3DplusT(StringContainer) - re-sort spatially sorted files (calls SortSeriesSlices unless explicitly told not to) - count how many images are at the position of the first image - assume this many 3D image blocks - for each block iterate all files in input, always skip "expected number of image blocks" - count how many images can be sorted into the first block - expect this many images for all other blocks - verify expectation, note image series as "cannot load as 3D+t" if expectation not met DataNode LoadDicom(StringContainer) - exception for Philips3D - if either ( - sorting by SortIntoBlocksFor3DplusT cannot form a valid 3D+t block OR - method is invoked to not load 3D+t ) - load only first block resulting from SortIntoBlocksFor3DplusT - else - create mitk::Image with multiple time step, assign result of loading each image block to one time step */ class MITK_CORE_EXPORT DicomSeriesReader { public: /** \brief For lists of filenames */ typedef std::vector StringContainer; typedef std::pair TwoStringContainers; /** \brief For multiple lists of filenames, assigned an ID each */ typedef std::map UidFileNamesMap; /** \brief Interface for the progress callback */ typedef void (*UpdateCallBackMethod)(float); /** \brief Checks if a specific file contains DICOM information. */ static bool IsDicom(const std::string &filename); /** \brief Find all series (and sub-series -- see details) in a particular directory. For each series, an enumeration of the files contained in it is created. \return The resulting map will map SeriesInstanceUID to UNSORTED lists of file names. SeriesInstanceUID will be CHANGED to be unique for each set of file names that should be 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. */ static UidFileNamesMap GetSeries(const std::string &dir, 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, const StringContainer &restrictions = StringContainer()); /** \brief see other GetSeries(). Instead of a directory name, provide an explicit list of file names. */ static UidFileNamesMap + GetSeries(const StringContainer& files, bool sortTo3DPlust, const StringContainer &restrictions = StringContainer()); + + /** + \brief Legacy interface. + + Use GetSeries(const StringContainer& files, bool sortTo3DPlust, const StringContainer &restrictions) instead. + */ + static + UidFileNamesMap GetSeries(const StringContainer& files, const StringContainer &restrictions = StringContainer()); /** relevant code that is matched here is in itkImageSeriesReader.txx (ImageSeriesReader::GenerateOutputInformation(void)), lines 176 to 245 (as of ITK 3.20) */ static TwoStringContainers AnalyzeFileForITKImageSeriesReaderSpacingAssumption(const StringContainer& files, const gdcm::Scanner::MappingType& tagValueMappings); /** 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); /** 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, 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, UpdateCallBackMethod callback = 0); /** \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(); protected: #if GDCM_MAJOR_VERSION >= 2 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 ); static std::string CreateSeriesIdentifierPart( gdcm::Scanner::TagToValue& tagValueMap, const gdcm::Tag& tag ); static std::string IDifyTagValue(const std::string& value); #endif #ifdef MITK_USE_GDCMIO typedef itk::GDCMImageIO DcmIoType; #else typedef itk::DICOMImageIO2 DcmIoType; #endif /** \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 Performs the 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, 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&, CallbackCommand* command = NULL); #if GDCM_MAJOR_VERSION >= 2 /** \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. */ static std::list SortIntoBlocksFor3DplusT( const StringContainer& presortedFilenames, 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); #endif }; } #endif /* MITKDICOMSERIESREADER_H_ */