diff --git a/Core/Code/IO/mitkDicomSeriesReader.cpp b/Core/Code/IO/mitkDicomSeriesReader.cpp index 6c7389be12..ace3c827a8 100644 --- a/Core/Code/IO/mitkDicomSeriesReader.cpp +++ b/Core/Code/IO/mitkDicomSeriesReader.cpp @@ -1,596 +1,893 @@ /*========================================================================= 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. =========================================================================*/ +// uncomment for learning more about the internal sorting mechanisms +//#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_ERROR << "Calling LoadDicomSeries with empty filename string container. Aborting."; - return false; + 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::UidFileNamesMap -DicomSeriesReader::GetSeries(const std::string &dir, const StringContainer &restrictions) +DicomSeriesReader::TwoStringContainers +DicomSeriesReader::AnalyzeFileForITKImageSeriesReaderSpacingAssumption( + const StringContainer& files, + const gdcm::Scanner::MappingType& tagValueMappings_) { - UidFileNamesMap map; // result variable + // result.first = files that fit ITK's assumption + // result.second = files that do not fit, should be run through AnalyzeFileForITKImageSeriesReaderSpacingAssumption() again + 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); + 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] << ")"; + + if ( lastOriginInitialized && (thisOrigin == lastOrigin) ) + { + MITK_DEBUG << "Sort away " << *fileIter << " for separate time step"; // we already have one occupying this position + result.second.push_back( *fileIter ); + } + 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 + + // 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 well, but it is just a name at this point + std::string thisOrientationString = tagValueMappings[fileIter->c_str()][tagImageOrientation]; + + 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); + + if (distance > 0.001) // mitk::eps is too small; 1/1000 of a mm should be enough to detect tilt + { + MITK_WARN << "Series seems to contain a tilted geometry. Will load series as many single slices."; + MITK_WARN << "Distance of expected slice origin from actual slice 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); + } + + // recored current origin for reference in later iterations + if ( !lastOriginInitialized || thisOrigin != lastOrigin ) + { + lastDifferentOrigin = thisOrigin; + } + + lastOrigin = thisOrigin; + lastOriginInitialized = true; + } + + 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) +{ /** 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 map; // preliminary result, refined into the final result mapOf3DPlusTBlocks #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 + // use GDCM directly, itk::GDCMSeriesFileNames does not work with GDCM 2 - // set directory - gdcm::Directory gDirectory; - gDirectory.Load( dir.c_str(), false ); // non-recursive - const gdcm::Directory::FilenamesType &gAllFiles = gDirectory.GetFilenames(); + // 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) - // scann for relevant tags in dicom files + // 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( gDirectory.GetFilenames() ) ) + if ( !scanner.Scan( files ) ) { - MITK_ERROR << "gdcm::Scanner failed scanning " << dir; + 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 << "File " << fileIter->first << std::endl; + 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 (i.e. blocks of equal z-spacing), + // 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 ) + { + map[ groupIter->first ] = SortSeriesSlices( groupIter->second ); // sort each slice group spatially + } + + UidFileNamesMap mapOf3DPlusTBlocks; // final result of this function + 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() ); + + // 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 + } + + // 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 + + // 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 + // TODO should be safe, but a little try/catch or other error handling wouldn't hurt + std::string thisOriginString = scanner.GetValue( mapOf3DBlocks[thisBlockKey].front().c_str(), tagImagePositionPatient ); + std::string previousOriginString = scanner.GetValue( mapOf3DBlocks[previousBlockKey].front().c_str(), tagImagePositionPatient ); + + // also compare last origin, because this might differ if z-spacing is different + std::string thisDestinationString = scanner.GetValue( mapOf3DBlocks[thisBlockKey].back().c_str(), tagImagePositionPatient ); + std::string previousDestinationString = scanner.GetValue( mapOf3DBlocks[previousBlockKey].back().c_str(), tagImagePositionPatient ); + + bool identicalOrigins( (thisOriginString == previousOriginString) && (thisDestinationString == previousDestinationString) ); + + 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 i = map.begin(); i != map.end(); ++i ) + for ( UidFileNamesMap::const_iterator groupIter = map.begin(); groupIter != map.end(); ++groupIter ) { - MITK_DEBUG << "Entry " << i->first << " with " << i->second.size() << " files"; + 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 '.' - MITK_DEBUG << "ID: " << constructedID; 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 61210dc06c..ed27355117 100644 --- a/Core/Code/IO/mitkDicomSeriesReader.h +++ b/Core/Code/IO/mitkDicomSeriesReader.h @@ -1,311 +1,441 @@ /*========================================================================= 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. + \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 + + \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 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. + 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: - After loading the series, the generated DataNode is not named. The caller is - responsible for given the node a name. + 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 directory known here: /home/me/dicom + // only a directory is known at this point: /home/who/dicom - DicomSeriesReader::UidFileNamesMap allImageBlocks = DicomSeriesReader::GetSeries("/home/me/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. - // optional step: sorting - std::string seriesToLoad = allImageBlocks[...]; // decide for yourself - - DicomSeriesReader::StringContainer oneBlockSorted = DicomSeriesReader::SortSeriesSlices( seriesToLoad ); + 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 ); // would sort again, but we could set the sort flag to false + DataNode::Pointer node = DicomSeriesReader::LoadDicomSeries( oneBlockSorted ); - Image::Pointer image = dynamic_cast( node.GetData() ); + 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 + - (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) - \section DicomSeriesReader_sorting Sorting into 3D+t blocks + \subsection DicomSeriesReader_sorting2 Step 2: Sort slices spatially - TODO: describe strategy of sorting images into blocks of 3D image stacks + 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) - \section DicomSeriesReader_limitations Assumptions and limitations + \subsection DicomSeriesReader_sorting3 Step 3: Ensure equal z spacing - 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. + 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. - \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") + 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. If an unexpected gap is detected, the block is split up. - \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) + Slices that share a position in space are also separated 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). */ class MITK_CORE_EXPORT DicomSeriesReader { public: /** - \brief For lists of filenames + \brief Lists of filenames. */ typedef std::vector StringContainer; /** - \brief For multiple lists of filenames, assigned an ID each + \brief For grouped lists of filenames, assigned an ID each. */ typedef std::map UidFileNamesMap; /** - \brief Interface for the progress callback + \brief Interface for the progress callback. */ typedef void (*UpdateCallBackMethod)(float); /** - \brief Checks if a specific file contains DICOM information. + \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 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. + \brief see other GetSeries(). + + Find all series (and sub-series -- see details) in a particular directory. */ 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()); /** - 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 + \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 StringContainer SortSeriesSlices(const StringContainer &unsortedFilenames); + static + UidFileNamesMap + GetSeries(const StringContainer& files, bool sortTo3DPlust, 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, 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, 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); +protected: + /** - \brief Provide combination of preprocessor defines that was active during compilation. + \brief for internal sorting. + */ + typedef std::pair TwoStringContainers; - 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. + /** + \brief Ensure an equal z-spacing for a group of files. + + 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 std::string GetConfigurationString(); + static + TwoStringContainers + AnalyzeFileForITKImageSeriesReaderSpacingAssumption(const StringContainer& files, const gdcm::Scanner::MappingType& tagValueMappings); + + /** + \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); -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 ); + + /** + \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); #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. + \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, 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. + + \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, 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_ */ diff --git a/Core/Code/Testing/CMakeLists.txt b/Core/Code/Testing/CMakeLists.txt index 2e3e2808a3..cbff306220 100644 --- a/Core/Code/Testing/CMakeLists.txt +++ b/Core/Code/Testing/CMakeLists.txt @@ -1,26 +1,29 @@ MITK_CREATE_MODULE_TESTS(LABELS MITK-Core) # MITK_INSTALL_TARGETS(EXECUTABLES MitkTestDriver) mitkAddCustomModuleTest(mitkPicFileReaderTest_emptyFile mitkPicFileReaderTest ${CMAKE_CURRENT_SOURCE_DIR}/Data/emptyFile.pic) mitkAddCustomModuleTest(mitkPicFileReaderTest_emptyGzipFile mitkPicFileReaderTest ${CMAKE_CURRENT_SOURCE_DIR}/Data/emptyFile.pic.gz) mitkAddCustomModuleTest(mitkDICOMLocaleTest_spacingOk_CT mitkDICOMLocaleTest ${MITK_DATA_DIR}/spacing-ok-ct.dcm) mitkAddCustomModuleTest(mitkDICOMLocaleTest_spacingOk_MR mitkDICOMLocaleTest ${MITK_DATA_DIR}/spacing-ok-mr.dcm) mitkAddCustomModuleTest(mitkDICOMLocaleTest_spacingOk_SC mitkDICOMLocaleTest ${MITK_DATA_DIR}/spacing-ok-sc.dcm) mitkAddCustomModuleTest(mitkEventMapperTest_Test1And2 mitkEventMapperTest ${MITK_DATA_DIR}/TestStateMachine1.xml ${MITK_DATA_DIR}/TestStateMachine2.xml) mitkAddCustomModuleTest(mitkNodeDependentPointSetInteractorTest mitkNodeDependentPointSetInteractorTest ${MITK_DATA_DIR}/Pic3D.pic.gz ${MITK_DATA_DIR}/BallBinary30x30x30.pic.gz) mitkAddCustomModuleTest(mitkDataStorageTest_US4DCyl mitkDataStorageTest ${MITK_DATA_DIR}/US4DCyl.pic.gz) mitkAddCustomModuleTest(mitkStateMachineFactoryTest_TestStateMachine1_2 mitkStateMachineFactoryTest ${MITK_DATA_DIR}/TestStateMachine1.xml ${MITK_DATA_DIR}/TestStateMachine2.xml) ADD_TEST(mitkPointSetLocaleTest ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/${TESTDRIVER} mitkPointSetLocaleTest ${MITK_DATA_DIR}/pointSet.mps) SET_PROPERTY(TEST mitkPointSetLocaleTest PROPERTY LABELS MITK-Core) ADD_TEST(mitkPicFileReaderTest_emptyGzipFile ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/${TESTDRIVER} mitkPicFileReaderTest ${CMAKE_CURRENT_SOURCE_DIR}/Data/emptyFile.pic.gz) SET_PROPERTY(TEST mitkPicFileReaderTest_emptyGzipFile PROPERTY LABELS MITK-Core) ADD_TEST(mitkImageTest_brainImage ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/${TESTDRIVER} mitkImageTest ${MITK_DATA_DIR}/brain.mhd) SET_PROPERTY(TEST mitkImageTest_brainImage PROPERTY LABELS MITK-Core) + +add_subdirectory(DICOMTesting) + diff --git a/Core/Code/Testing/DICOMTesting/CMakeLists.txt b/Core/Code/Testing/DICOMTesting/CMakeLists.txt new file mode 100644 index 0000000000..2ba6f34a77 --- /dev/null +++ b/Core/Code/Testing/DICOMTesting/CMakeLists.txt @@ -0,0 +1,35 @@ +if (BUILD_TESTING) +if (MITK_USE_GDCMIO) +if (GDCM_DIR) + +# clear variables from prior files.cmake +set(MODULE_TESTS) +set(MODULE_IMAGE_TESTS) +set(MODULE_TESTIMAGES) +set(MODULE_CUSTOM_TESTS) +set(H_FILES) +set(CPP_FILES) + +# now create a new module only for testing purposes +MITK_CREATE_MODULE( mitkDICOMTesting + DEPENDS Mitk # mitkCore.so +) + + +# add helpful applications +MITK_USE_MODULE( mitkDICOMTesting ) +include_directories(${CMAKE_CURRENT_SOURCE_DIR} ${CMAKE_CURRENT_BINARY_DIR} ${ALL_INCLUDE_DIRECTORIES}) + +# dumps out image information +add_executable(DumpDICOMMitkImage DumpDICOMMitkImage.cpp) +target_link_libraries(DumpDICOMMitkImage ${ALL_LIBRARIES}) + +# compares dumped out image information against reference dump +add_executable(VerifyDICOMMitkImageDump VerifyDICOMMitkImageDump.cpp) +target_link_libraries(VerifyDICOMMitkImageDump ${ALL_LIBRARIES}) + +add_subdirectory(Testing) + +endif() +endif() +endif() diff --git a/Core/Code/Testing/DICOMTesting/DICOMTesting.dox b/Core/Code/Testing/DICOMTesting/DICOMTesting.dox new file mode 100644 index 0000000000..a751217c9c --- /dev/null +++ b/Core/Code/Testing/DICOMTesting/DICOMTesting.dox @@ -0,0 +1,144 @@ +/** + +\page DICOMTesting MITK DICOM testing + +\section DICOMTesting_introduction Introduction + +Reading DICOM data into mitk::Images is a complicated business since DICOM and MITK have very different ideas of images. + +DicomSeriesReader brings DICOM and MITK as close together as possible by offering methods to load DICOM CT and MR images +into mitk::Images, optionally grouping slices to 3D+t images. + +Since there are so many possible sources for mistakes with any change to this loading process, testing the many +assumptions implemented in DicomSeriesReader is worthwhile. This document describes what kind of tests are implemented +and what how. + +\section DICOMTesting_problem Problem description + +The task of loading DICOM files into mitk::Images is a challenge because of major differences +in the way that DICOM and MITK represent images: + + - DICOM images + - are mostly stored as one slice per file + - do not describe how they can be arraged into image stacks with orthogonal axis + - sometimes they cannot be arranged in image stacks as described above (e.g. tilted gantry) + - mitk::Image (at least its mature areas) + - represents image stacks with orthogonal axis (nothing like a tilted gantry) + - have a concept of a fourth dimension (time steps) + +Because image processing applications based on MITK still need a way to load DICOM images into +mitk::Image objects, MITK needs to put a lot of effort into building image stacks as best as it can. And this needs to be well tested. + +For more background information, see David Clunie's most valuable posts on comp.protocols.dicom, e.g.: + + - http://groups.google.com/group/comp.protocols.dicom/browse_thread/thread/6db91972e161f0d4/6e0304ac264a6eb5 + - http://groups.google.com/group/comp.protocols.dicom/browse_thread/thread/e9bd1497bea3e66b/187a7dc8810613d2 + - http://groups.google.com/group/comp.protocols.dicom/browse_thread/thread/5d80bb0b7fafcb81/cf96119e3b024ed8 + - http://groups.google.com/group/comp.protocols.dicom/browse_thread/thread/4568635e083a3fba/e2a8ceec23032601 + +\section DICOMTesting_testidea Test principle + +The general idea for DICOM loaing tests is to run a set of known DICOM files through DicomSeriesReader's methods +GetSeries() and LoadDicomSeries() to generate mitk::Images. These images are then compared to expected image properties, +such as the number of individual mitk::Images, positions, orientations, spacings, etc. + +Stored expectations look like this (should be self-explanatory): +\verbatim +-- Image 1 +Pixeltype: s +BitsPerPixel: 16 +Dimension: 4 +Dimensions: 64 64 6 1 +Geometry: + Matrix: 5.25 0 0 0 5.2468 0.139598 0 -0.183222 3.99757 + Offset: -159.672 -309.974 -69.0122 + Center: 0 0 0 + Translation: -159.672 -309.974 -69.0122 + Scale: 1 1 1 + Origin: -159.672 -309.974 -69.0122 + Spacing: 5.25 5.25 4 + TimeBounds: 0 1 + +-- Image 2 +Pixeltype: s +BitsPerPixel: 16 +Dimension: 4 +Dimensions: 64 64 41 1 +Geometry: + Matrix: 5.25 0 0 0 5.25 0 0 0 4 + Offset: -160.672 -311.672 -285 + Center: 0 0 0 + Translation: -160.672 -311.672 -285 + Scale: 1 1 1 + Origin: -160.672 -311.672 -285 + Spacing: 5.25 5.25 4 + TimeBounds: 0 1 +\endverbatim + + + +\section DICOMTesting_implementation Implementation + +\section DICOMTesting_implementation_utils Test helpers (applications and classes) + +Application DumpDICOMMitkImage + +Takes a list of DICOM images, loads them using TestDICOMLoading, then dumps information +about the resulting mitk::Images to standard output. + +This application is helpful when defining reference data for tests. + +Application VerifyDICOMMitkImageDump + +Takes a list of DICOM images and loads them using TestDICOMLoading. +Takes a dump file as generated by DumpDICOMMitkImage, parses it and +compares it to actually generated mitk::Images. + +This application is used to implement the majority of test cases. They all +load images, then verify the expected result structure. + +Class TestDICOMLoading + +\section PageDICOMLoadingTests_testcaseimplementation Test case implementation + +Individual test cases are stored in the MITK-Data repository and constructed by Core/Code/Testing/DICOMTesting/Testing/CMakeLists.txt + +The CMake file parses given directories for subdirectories containing specific test cases. Each directory contains two files: + - File "input": lists DICOM files that should be loaded for this test case + - File "expected.dump": contains the image properties in the above mentioned dump format + +Each test case is translated into a CTest test which evaluates the return value of a call to VerifyDICOMMitkImageDump. + +\section PageDICOMLoadingTests_testcases Implemented test cases + +From test set TinyCTAbdomen (see description.txt in this directory for details on test images): + + - singleslice : just a single slice (should work and contain meaningful spacing) + - two_slices : two slices, spacing should be calculated correctly + - all : load a "normal" series as a single 3D block + - 3D_and_T : load a small set of slices with multiple time-steps + - diff_orientation : load a set of files containing two differently oriented image blocks + - diff_orientation_gaps : load a set of files containing two differently oriented image blocks, each missing slices, so blocks must be split + - diff_spacing : load a set of files containint two set of slices with different spacings + - gap : load slices that cannot form a single 3D block, because single slices are missing + - gaps : slices missing in differnt locations, so multiple splits needed + - unsorted_gaps : as before, just file names are unsorted + - single_negative_spacing : from reported bug related to single MR images with misleading spacing information + - tilted_gantry : slice origins do not align along first slice normal (happens with tilted gantries) + + +\section DICOMTesting_othertests Specific other tests + +This list is meant to provide an up-to-date list of all implemented DICOM loading tests. +If you ever find this outdated, please update it or make the persons who invalidated the list update it. + +mitkDICOMTestingSanityTest_* +These tests implement basic testing of the implemented helper classes. The tests use DicomSeriesReader to load +a number of DICOM image. They verify: + - DicomSeriesReader recognizes all input files as DICOM images + - DicomSeriesReader generates a number of mitk::Images from the DICOM images + - the number of mitk::Images matches a number given on the command line or CTest's add_test() + - helper methods in class TestDICOMLoading make minimal sense (comparison of an image dump to itself must be valid) + +*/ + diff --git a/Core/Code/Testing/DICOMTesting/DumpDICOMMitkImage.cpp b/Core/Code/Testing/DICOMTesting/DumpDICOMMitkImage.cpp new file mode 100644 index 0000000000..a6f0d63b2b --- /dev/null +++ b/Core/Code/Testing/DICOMTesting/DumpDICOMMitkImage.cpp @@ -0,0 +1,23 @@ +#include "mitkTestDICOMLoading.h" + +int main(int argc, char** argv) +{ + mitk::TestDICOMLoading loader; + mitk::TestDICOMLoading::StringContainer files; + + for (int arg = 1; arg < argc; ++arg) files.push_back( argv[arg] ); + + mitk::TestDICOMLoading::ImageList images = loader.LoadFiles(files); + + // combine individual dumps in a way that VerifyDICOMMitkImageDump is able to separate again. + // I.e.: when changing this piece of code, always change VerifyDICOMMitkImageDump, too. + unsigned int imageCounter(0); + for ( mitk::TestDICOMLoading::ImageList::const_iterator imageIter = images.begin(); + imageIter != images.end(); + ++imageIter ) + { + std::cout << "-- Image " << ++imageCounter << "\n"; + std::cout << loader.DumpImageInformation( *imageIter ) << "\n"; + } +} + diff --git a/Core/Code/Testing/DICOMTesting/Testing/CMakeLists.txt b/Core/Code/Testing/DICOMTesting/Testing/CMakeLists.txt new file mode 100644 index 0000000000..72de94f8bc --- /dev/null +++ b/Core/Code/Testing/DICOMTesting/Testing/CMakeLists.txt @@ -0,0 +1,59 @@ + +MITK_CREATE_MODULE_TESTS(LABELS MITK-Core) + + +include(mitkFunctionAddTestLabel) + +# verify minimum expectations: +# files are recognized as DICOM +# loading files results in a given number of images +mitkAddCustomModuleTest(mitkDICOMTestingSanityTest_NoFiles mitkDICOMTestingSanityTest 0) +mitkAddCustomModuleTest(mitkDICOMTestingSanityTest_CTImage mitkDICOMTestingSanityTest 1 ${MITK_DATA_DIR}/spacing-ok-ct.dcm) +mitkAddCustomModuleTest(mitkDICOMTestingSanityTest_MRImage mitkDICOMTestingSanityTest 1 ${MITK_DATA_DIR}/spacing-ok-mr.dcm) +mitkAddCustomModuleTest(mitkDICOMTestingSanityTest_SCImage mitkDICOMTestingSanityTest 1 ${MITK_DATA_DIR}/spacing-ok-sc.dcm) + +set(VERIFY_DUMP_CMD ${CMAKE_RUNTIME_OUTPUT_DIRECTORY}/VerifyDICOMMitkImageDump) + +set(CT_ABDOMEN_DIR ${MITK_DATA_DIR}/TinyCTAbdomen) + +# these variables could be passed as parameters to a generic test creation function +set(CURRENT_DATASET_DIR ${CT_ABDOMEN_DIR}) +set(TESTS_DIR Tests) +set(INPUT_LISTNAME input) +set(EXPECTED_DUMP expected.dump) + +file(GLOB_RECURSE allInputs ${CURRENT_DATASET_DIR}/${TESTS_DIR}/*/${INPUT_LISTNAME}) + +function(expectFileExists filename) + if (NOT EXISTS ${filename}) + message(SEND_ERROR "Test case expected file ${filename} which does not exist! Please fix your CMake code or file layout.") + endif (NOT EXISTS ${filename}) +endfunction(expectFileExists) + +foreach(input ${allInputs}) + string(REGEX REPLACE ".*${TESTS_DIR}/([^/]+)/.*" "\\1" input ${input}) + set(inputfilelist "${CURRENT_DATASET_DIR}/${TESTS_DIR}/${input}/${INPUT_LISTNAME}") + set(expecteddump "${CURRENT_DATASET_DIR}/${TESTS_DIR}/${input}/${EXPECTED_DUMP}") + set(testname "DICOM_Load_${input}") + + message(STATUS "DICOM loading test case '${input}'") + expectFileExists(${inputfilelist}) + expectFileExists(${expecteddump}) + + # TODO provide error messages if input not valid + + set(dicomFiles) # clear list + file(STRINGS ${inputfilelist} rawDicomFiles) + foreach(raw ${rawDicomFiles}) + set(fileWithFullPath ${CURRENT_DATASET_DIR}/${raw}) + list(APPEND dicomFiles ${fileWithFullPath}) + endforeach(raw ${rawDicomFiles}) + + #message(STATUS " Load ${dicomFiles}") + add_test(${testname} + ${VERIFY_DUMP_CMD} ${expecteddump} ${dicomFiles}) + mitkFunctionAddTestLabel(${testname}) + +endforeach(input allInputs) + + diff --git a/Core/Code/Testing/DICOMTesting/Testing/files.cmake b/Core/Code/Testing/DICOMTesting/Testing/files.cmake new file mode 100644 index 0000000000..92af05e67b --- /dev/null +++ b/Core/Code/Testing/DICOMTesting/Testing/files.cmake @@ -0,0 +1,6 @@ + +# tests with no extra command line parameter +SET(MODULE_CUSTOM_TESTS + mitkDICOMTestingSanityTest.cpp +) + diff --git a/Core/Code/Testing/DICOMTesting/Testing/mitkDICOMTestingSanityTest.cpp b/Core/Code/Testing/DICOMTesting/Testing/mitkDICOMTestingSanityTest.cpp new file mode 100644 index 0000000000..08262dc82b --- /dev/null +++ b/Core/Code/Testing/DICOMTesting/Testing/mitkDICOMTestingSanityTest.cpp @@ -0,0 +1,54 @@ + + +#include "mitkTestDICOMLoading.h" +#include "mitkTestingMacros.h" + +int mitkDICOMTestingSanityTest(int argc, char** const argv) +{ + MITK_TEST_BEGIN("DICOMTestingSanity") + + mitk::TestDICOMLoading loader; + mitk::TestDICOMLoading::StringContainer files; + + // adapt expectations depending on configuration + std::string configuration = mitk::DicomSeriesReader::GetConfigurationString(); + MITK_TEST_OUTPUT(<< "Configuration: " << configuration) + /* + MITK_TEST_CONDITION_REQUIRED( configuration.find( "GDCM_VERSION: 2." ) != std::string::npos, + "Expect at least GDCM version 2" ) + */ + + // load files from commandline + int numberOfExpectedImages = 0; + if (argc > 1) numberOfExpectedImages = atoi(argv[1]); + for (int arg = 2; arg < argc; ++arg) files.push_back( argv[arg] ); + + // verify all files are DICOM + for (mitk::TestDICOMLoading::StringContainer::const_iterator fileIter = files.begin(); + fileIter != files.end(); + ++fileIter) + { + MITK_TEST_CONDITION_REQUIRED( mitk::DicomSeriesReader::IsDicom(*fileIter) , *fileIter << " is recognized as loadable DICOM object" ) + + } + + // compare with expected number of images from commandline + mitk::TestDICOMLoading::ImageList images = loader.LoadFiles(files); + MITK_TEST_CONDITION_REQUIRED( images.size() == numberOfExpectedImages, "Loading " << files.size() + << " files from commandline results in " << numberOfExpectedImages + << " images (see test invocation)" ) + + // check dump equality (dumping image information must always equal itself) + for ( mitk::TestDICOMLoading::ImageList::const_iterator imageIter = images.begin(); + imageIter != images.end(); + ++imageIter ) + { + const mitk::Image* image = *imageIter; + MITK_TEST_CONDITION( loader.CompareImageInformationDumps( loader.DumpImageInformation(image), + loader.DumpImageInformation(image) ) == true, + "Image information dumping is able to reproduce its result." ) + } + + MITK_TEST_END() +} + diff --git a/Core/Code/Testing/DICOMTesting/VerifyDICOMMitkImageDump.cpp b/Core/Code/Testing/DICOMTesting/VerifyDICOMMitkImageDump.cpp new file mode 100644 index 0000000000..2bd6ee1746 --- /dev/null +++ b/Core/Code/Testing/DICOMTesting/VerifyDICOMMitkImageDump.cpp @@ -0,0 +1,107 @@ + +#include "mitkTestDICOMLoading.h" + +std::vector LoadDumps(const std::string& fileName) +{ + std::vector separatedDumps; + + std::ifstream fileStream( fileName.c_str() ); + + std::string buffer; + std::string line; + while(fileStream){ + std::getline(fileStream, line); + + if (line.find("-- Image ") == 0) + { + // separator: starts a new image block + if ( !buffer.empty() ) + { + // unless this is the first block + separatedDumps.push_back(buffer); + buffer.clear(); + } + } + else + { + buffer += line + "\n"; + } + } + fileStream.close(); + + // eat last image dump + if ( !buffer.empty() ) + { + separatedDumps.push_back(buffer); + buffer.clear(); + } + + return separatedDumps; +} + + +int main(int argc, char** argv) +{ + /** + Loads a list of DICOM images, compares generated mitk::Images against stored references. + + first argument: file with reference dumps + following arguments: file names to load + */ + + if (argc < 2) + { + MITK_ERROR << "Usage:"; + MITK_ERROR << " " << argv[0] << " reference.dump file1 [file2 .. fileN]"; + MITK_ERROR << " "; + MITK_ERROR << " Loads all DICOM images in file1 to fileN as MITK images "; + MITK_ERROR << " and compares loaded images against stored expectations (dumps)."; + MITK_ERROR << " See also DumpDICOMMitkImage (generates dumps)"; + return EXIT_FAILURE; + } + + mitk::TestDICOMLoading loader; + mitk::TestDICOMLoading::StringContainer files; + + // TODO load reference dumps + + // load from argv[1] + // separate at "-- Image n" + // store in an array of dumps + + std::vector expectedDumps = LoadDumps(argv[1]); + + for (int arg = 2; arg < argc; ++arg) files.push_back( argv[arg] ); + + mitk::TestDICOMLoading::ImageList images = loader.LoadFiles(files); + + unsigned int imageCounter(0); + for ( mitk::TestDICOMLoading::ImageList::const_iterator imageIter = images.begin(); + imageIter != images.end(); + ++imageIter, ++imageCounter ) + { + std::string imageDump = loader.DumpImageInformation( *imageIter ); + + if (imageCounter >= expectedDumps.size()) + { + MITK_ERROR << "Loader produces more images than expected. Aborting after image " << (imageCounter-1); + MITK_INFO << "Image " << imageCounter << " loaded as:\n" << imageDump; + return EXIT_FAILURE; + } + + bool loadedAsExpected = loader.CompareImageInformationDumps( expectedDumps[imageCounter], imageDump ); + if (loadedAsExpected) + { + MITK_INFO << "Image " << imageCounter << " loads as expected."; + } + else + { + MITK_ERROR << "Image " << imageCounter << " did not load as expected."; + MITK_INFO << "Expected: \n" << expectedDumps[imageCounter] << "\nGot:\n" << imageDump; + return EXIT_FAILURE; + } + } + + return EXIT_SUCCESS; +} + diff --git a/Core/Code/Testing/DICOMTesting/files.cmake b/Core/Code/Testing/DICOMTesting/files.cmake new file mode 100644 index 0000000000..f6fd6ff9d6 --- /dev/null +++ b/Core/Code/Testing/DICOMTesting/files.cmake @@ -0,0 +1,4 @@ + +SET(CPP_FILES + mitkTestDICOMLoading.cpp +) diff --git a/Core/Code/Testing/DICOMTesting/mitkTestDICOMLoading.cpp b/Core/Code/Testing/DICOMTesting/mitkTestDICOMLoading.cpp new file mode 100644 index 0000000000..3319055d7f --- /dev/null +++ b/Core/Code/Testing/DICOMTesting/mitkTestDICOMLoading.cpp @@ -0,0 +1,331 @@ +//#define MBILOG_ENABLE_DEBUG + +#include "mitkTestDICOMLoading.h" + +#include + +mitk::TestDICOMLoading::TestDICOMLoading() +:m_PreviousCLocale(NULL) +{ +} + +void mitk::TestDICOMLoading::SetDefaultLocale() +{ + // remember old locale only once + if (m_PreviousCLocale == NULL) + { + m_PreviousCLocale = setlocale(LC_NUMERIC, NULL); + + // set to "C" + setlocale(LC_NUMERIC, "C"); + + m_PreviousCppLocale = std::cin.getloc(); + + std::locale l( "C" ); + std::cin.imbue(l); + std::cout.imbue(l); + } +} + +void mitk::TestDICOMLoading::ResetUserLocale() +{ + if (m_PreviousCLocale) + { + setlocale(LC_NUMERIC, m_PreviousCLocale); + + std::cin.imbue(m_PreviousCppLocale); + std::cout.imbue(m_PreviousCppLocale); + + m_PreviousCLocale = NULL; + } +} + + + +mitk::TestDICOMLoading::ImageList mitk::TestDICOMLoading::LoadFiles( const StringContainer& files ) +{ + for (StringContainer::const_iterator iter = files.begin(); + iter != files.end(); + ++iter) + { + MITK_DEBUG << "File " << *iter; + } + + ImageList result; + + DicomSeriesReader::UidFileNamesMap seriesInFiles = DicomSeriesReader::GetSeries( files ); + + // TODO sort series UIDs, implementation of map iterator might differ on different platforms (or verify this is a standard topic??) + for (DicomSeriesReader::UidFileNamesMap::const_iterator seriesIter = seriesInFiles.begin(); + seriesIter != seriesInFiles.end(); + ++seriesIter) + { + StringContainer files = seriesIter->second; + + DataNode::Pointer node = DicomSeriesReader::LoadDicomSeries( files ); + + if (node.IsNotNull()) + { + Image::Pointer image = dynamic_cast( node->GetData() ); + + result.push_back( image ); + } + else + { + } + } + + return result; +} + +std::string +mitk::TestDICOMLoading::TypeIDToString(const std::type_info& ti) +{ + if (ti == typeid(unsigned char)) + return "UCHAR"; + else if (ti == typeid(char)) + return "CHAR"; + else if (ti == typeid(unsigned short)) + return "USHORT"; + else if (ti == typeid(short)) + return "SHORT"; + else if (ti == typeid(unsigned int)) + return "UINT"; + else if (ti == typeid(int)) + return "INT"; + else if (ti == typeid(long unsigned int)) + return "ULONG"; + else if (ti == typeid(long int)) + return "LONG"; + else if (ti == typeid(float)) + return "FLOAT"; + else if (ti == typeid(double)) + return "DOUBLE"; + else + return "UNKNOWN"; +} + + +// add a line to stringstream result (see DumpImageInformation +#define DumpLine(field, data) DumpILine(0, field, data) + +// add an indented(!) line to stringstream result (see DumpImageInformation +#define DumpILine(indent, field, data) \ +{ \ + std::string DumpLine_INDENT; DumpLine_INDENT.resize(indent, ' ' ); \ + result << DumpLine_INDENT << field << ": " << data << "\n"; \ +} + +std::string +mitk::TestDICOMLoading::DumpImageInformation( const Image* image ) +{ + std::stringstream result; + + if (image == NULL) return result.str(); + + SetDefaultLocale(); + + // basic image data + DumpLine( "Pixeltype", TypeIDToString( *(image->GetPixelType().GetTypeId()) )); + DumpLine( "BitsPerPixel", image->GetPixelType().GetBpe() ); + DumpLine( "Dimension", image->GetDimension() ); + + result << "Dimensions: "; + for (unsigned int dim = 0; dim < image->GetDimension(); ++dim) + result << image->GetDimension(dim) << " "; + result << "\n"; + + // geometry data + result << "Geometry: \n"; + Geometry3D* geometry = image->GetGeometry(); + if (geometry) + { + AffineTransform3D* transform = geometry->GetIndexToWorldTransform(); + if (transform) + { + result << " " << "Matrix: "; + const AffineTransform3D::MatrixType& matrix = transform->GetMatrix(); + for (unsigned int i = 0; i < 3; ++i) + for (unsigned int j = 0; j < 3; ++j) + result << matrix[i][j] << " "; + result << "\n"; + + result << " " << "Offset: "; + const AffineTransform3D::OutputVectorType& offset = transform->GetOffset(); + for (unsigned int i = 0; i < 3; ++i) + result << offset[i] << " "; + result << "\n"; + + result << " " << "Center: "; + const AffineTransform3D::InputPointType& center = transform->GetCenter(); + for (unsigned int i = 0; i < 3; ++i) + result << center[i] << " "; + result << "\n"; + + result << " " << "Translation: "; + const AffineTransform3D::OutputVectorType& translation = transform->GetTranslation(); + for (unsigned int i = 0; i < 3; ++i) + result << translation[i] << " "; + result << "\n"; + + result << " " << "Scale: "; + const double* scale = transform->GetScale(); + for (unsigned int i = 0; i < 3; ++i) + result << scale[i] << " "; + result << "\n"; + + result << " " << "Origin: "; + const Point3D& origin = geometry->GetOrigin(); + for (unsigned int i = 0; i < 3; ++i) + result << origin[i] << " "; + result << "\n"; + + result << " " << "Spacing: "; + const Vector3D& spacing = geometry->GetSpacing(); + for (unsigned int i = 0; i < 3; ++i) + result << spacing[i] << " "; + result << "\n"; + + result << " " << "TimeBounds: "; + const TimeBounds timeBounds = geometry->GetTimeBounds(); + for (unsigned int i = 0; i < 2; ++i) + result << timeBounds[i] << " "; + result << "\n"; + + + } + } + + ResetUserLocale(); + + return result.str(); +} + +bool +mitk::TestDICOMLoading::CompareSpacedValueFields( const std::string& reference, + const std::string& test, + double eps ) +{ + // TODO this should better do a real comparison and tolerate float differences up to 'eps' + + return reference == test; +} + +bool +mitk::TestDICOMLoading::CompareImageInformationDumps( const std::string& referenceDump, + const std::string& testDump ) +{ + KeyValueMap reference = ParseDump(referenceDump); + KeyValueMap test = ParseDump(testDump); + + bool testResult(true); + + // verify all expected values + for (KeyValueMap::const_iterator refIter = reference.begin(); + refIter != reference.end(); + ++refIter) + { + const std::string& refKey = refIter->first; + const std::string& refValue = refIter->second; + + if ( test.find(refKey) != test.end() ) + { + const std::string& testValue = test[refKey]; + + MITK_DEBUG << refKey << ": " << refValue << " == " << testValue << "?"; + + testResult &= CompareSpacedValueFields( refValue, testValue ); + } + else + { + MITK_ERROR << "Reference dump contains a key'" << refKey << "' (value '" << refValue << "')." ; + MITK_ERROR << "This key is expected to be generated for tests (but was not). Most probably you need to update your test data."; + return false; + } + } + + // now check test dump does not contain any additional keys + for (KeyValueMap::const_iterator testIter = test.begin(); + testIter != test.end(); + ++testIter) + { + const std::string& key = testIter->first; + const std::string& value = testIter->second; + + if ( reference.find(key) == reference.end() ) + { + MITK_ERROR << "Test dump contains an unexpected key'" << key << "' (value '" << value << "')." ; + MITK_ERROR << "This key is not expected. Most probably you need to update your test data."; + return false; + } + } + + return testResult; +} + +mitk::TestDICOMLoading::KeyValueMap +mitk::TestDICOMLoading::ParseDump( const std::string& dump ) +{ + KeyValueMap parsedResult; + + std::string shredder(dump); + + std::stack surroundingKeys; + + std::stack expectedIndents; + expectedIndents.push(0); + + while (true) + { + std::string::size_type newLinePos = shredder.find( '\n' ); + if (newLinePos == std::string::npos || newLinePos == 0) break; + + std::string line = shredder.substr( 0, newLinePos ); + shredder = shredder.erase( 0, newLinePos+1 ); + + std::string::size_type keyPosition = line.find_first_not_of( ' ' ); + std::string::size_type colonPosition = line.find( ':' ); + + std::string key = line.substr(keyPosition, colonPosition - keyPosition); + std::string::size_type firstSpacePosition = key.find_first_of(" "); + if (firstSpacePosition != std::string::npos) + { + key.erase(firstSpacePosition); + } + + if ( keyPosition > expectedIndents.top() ) + { + // more indent than before + expectedIndents.push(keyPosition); + } + else if (keyPosition == expectedIndents.top() ) + { + if (!surroundingKeys.empty()) + { + surroundingKeys.pop(); // last of same length + } + } + else + { + // less indent than before + do expectedIndents.pop(); + while (expectedIndents.top() != keyPosition); // unwind until current indent is found + } + + if (!surroundingKeys.empty()) + { + key = surroundingKeys.top() + "." + key; // construct current key name + } + + surroundingKeys.push(key); // this is the new embracing key + + std::string value = line.substr(colonPosition+1); + + MITK_DEBUG << " Key: '" << key << "' value '" << value << "'" ; + + parsedResult[key] = value; // store parsing result + } + + return parsedResult; +} + diff --git a/Core/Code/Testing/DICOMTesting/mitkTestDICOMLoading.h b/Core/Code/Testing/DICOMTesting/mitkTestDICOMLoading.h new file mode 100644 index 0000000000..b261dc7a8f --- /dev/null +++ b/Core/Code/Testing/DICOMTesting/mitkTestDICOMLoading.h @@ -0,0 +1,64 @@ +#ifndef mitkTestDICOMLoading_h +#define mitkTestDICOMLoading_h + +#include "mitkDicomSeriesReader.h" + +#include "mitkDICOMTestingExports.h" + +namespace mitk +{ + +class mitkDICOMTesting_EXPORT TestDICOMLoading +{ + public: + + typedef DicomSeriesReader::StringContainer StringContainer; + typedef std::list NodeList; + typedef std::list ImageList; + + TestDICOMLoading(); + + ImageList + LoadFiles( const StringContainer& files ); + + /** + \brief Dump relevant image information for later comparison. + \sa CompareImageInformationDumps + */ + std::string + DumpImageInformation( const Image* image ); + + /** + \brief Compare two image information dumps. + \return true, if dumps are sufficiently equal (see parameters) + \sa DumpImageInformation + */ + bool + CompareImageInformationDumps( const std::string& reference, + const std::string& test ); + + private: + + typedef std::map KeyValueMap; + + void SetDefaultLocale(); + + void ResetUserLocale(); + + std::string TypeIDToString( const std::type_info& ); + + KeyValueMap ParseDump( const std::string& dump ); + + bool CompareSpacedValueFields( const std::string& reference, + const std::string& test, + double eps = mitk::eps ); + + const char* m_PreviousCLocale; + std::locale m_PreviousCppLocale; + +}; + +} + +#endif +