From de83970b84cb6117feb0904c7dd96ea16f755e76 Mon Sep 17 00:00:00 2001 From: tbaudier Date: Wed, 17 May 2017 16:54:57 +0200 Subject: [PATCH] Update GateSimulation2Dicom tool with test functions --- .../clitkGateSimulation2DicomGenericFilter.h | 6 + ...clitkGateSimulation2DicomGenericFilter.txx | 197 ++++++++++++++---- 2 files changed, 164 insertions(+), 39 deletions(-) diff --git a/tools/clitkGateSimulation2DicomGenericFilter.h b/tools/clitkGateSimulation2DicomGenericFilter.h index e3b78fe..7dbc725 100644 --- a/tools/clitkGateSimulation2DicomGenericFilter.h +++ b/tools/clitkGateSimulation2DicomGenericFilter.h @@ -110,8 +110,14 @@ namespace clitk bool m_Verbose; std::string m_InputFileName; + }; +//Copy dicom dictionary +void CopyDictionary (itk::MetaDataDictionary &fromDict, itk::MetaDataDictionary &toDict); + +//convert to std::string +template std::string NumberToString ( T Number ); } // end namespace clitk diff --git a/tools/clitkGateSimulation2DicomGenericFilter.txx b/tools/clitkGateSimulation2DicomGenericFilter.txx index f49bebc..39279dd 100644 --- a/tools/clitkGateSimulation2DicomGenericFilter.txx +++ b/tools/clitkGateSimulation2DicomGenericFilter.txx @@ -27,15 +27,21 @@ * ===================================================*/ +#include // clitk #include "clitkResampleImageWithOptionsFilter.h" #if GDCM_MAJOR_VERSION >= 2 #include "gdcmUIDGenerator.h" +#include +#include +#include #else #include "gdcmFile.h" #include "gdcmUtil.h" #endif +#include "itkImageRegionIterator.h" + namespace clitk { @@ -125,39 +131,37 @@ void GateSimulation2DicomGenericFilter::UpdateWithDimAndPixelType() { +#if GDCM_MAJOR_VERSION == 2 // ImageTypes typedef itk::Image InputImageType; typedef itk::Image OutputImageType; - - // Read the dicom directory typedef itk::ImageFileReader< InputImageType > ReaderType; + typedef itk::ImageSeriesReader< InputImageType > ReaderSeriesType; typedef itk::GDCMImageIO ImageIOType; - typedef itk::GDCMSeriesFileNames NamesGeneratorType; - - ImageIOType::Pointer gdcmIO = ImageIOType::New(); - std::string filename_out = m_ArgsInfo.outputFilename_arg; - // Output the dicom files - unsigned int numberOfFilenames = 1; - if (m_Verbose) { - std::cout << numberOfFilenames <<" were read in the directory "<= 2 + typename ReaderSeriesType::Pointer readerSeries = ReaderSeriesType::New(); + ImageIOType::Pointer gdcmIO = ImageIOType::New(); + std::string filename_out = m_ArgsInfo.outputFilename_arg; gdcmIO->LoadPrivateTagsOn(); gdcmIO->KeepOriginalUIDOn(); -#endif reader->SetImageIO( gdcmIO ); reader->SetFileName( m_ArgsInfo.inputModelFilename_arg ); + typename ReaderSeriesType::FileNamesContainer fileNames; + fileNames.push_back(m_ArgsInfo.inputModelFilename_arg); + readerSeries->SetImageIO( gdcmIO ); + readerSeries->SetFileNames( fileNames ); try { reader->Update(); + readerSeries->Update(); } catch (itk::ExceptionObject &excp) { - std::cerr << "Error: Exception thrown while reading the DICOM series!!" << std::endl; + std::cerr << "Error: Exception thrown while reading the DICOM model file !!" << std::endl; std::cerr << excp << std::endl; } + // Read the input (MHD file) typedef typename InputImageType::RegionType RegionType; typedef typename RegionType::SizeType SizeType; @@ -165,9 +169,11 @@ GateSimulation2DicomGenericFilter::UpdateWithDimAndPixelType() typename InputReaderType::Pointer volumeReader = InputReaderType::New(); volumeReader->SetFileName( m_InputFileName); volumeReader->Update(); - typename InputImageType::Pointer input = volumeReader->GetOutput(); - if (input->GetLargestPossibleRegion().GetSize() != reader->GetOutput()->GetLargestPossibleRegion().GetSize()) { + + +//TODO if the input size/spacing and dicom model ones are different +/* if (input->GetLargestPossibleRegion().GetSize() != reader->GetOutput()->GetLargestPossibleRegion().GetSize()) { // resampling is carried out on the fly if resolution or size between // the input mhd and input dicom series is different @@ -199,49 +205,162 @@ GateSimulation2DicomGenericFilter::UpdateWithDimAndPixelType() std::cerr << excp << std::endl; } } - -#if GDCM_MAJOR_VERSION < 2 - gdcmIO->SetKeepOriginalUID(true); -#endif +*/ //Read the dicom file to find the energy, the head & the steps (TODO: do it on the mhd filetext) + gdcm::Reader hreader; + hreader.SetFileName(m_ArgsInfo.inputModelFilename_arg); + hreader.Read(); + gdcm::DataSet& ds = hreader.GetFile().GetDataSet(); + int energyNumber, headNumber, rotationNumber; //Read the number of energy, of head and rotation steps + gdcm::Attribute<0x54,0x11> energyNumber_att; + energyNumber_att.SetFromDataSet(hreader.GetFile().GetDataSet()); + energyNumber = energyNumber_att.GetValue(); + + gdcm::Attribute<0x54,0x21> headNumber_att; + headNumber_att.SetFromDataSet(hreader.GetFile().GetDataSet()); + headNumber = headNumber_att.GetValue(); + + gdcm::Tag squenceTag(0x54,0x52); //rotation, read a sequence first + const gdcm::DataElement &sequence = hreader.GetFile().GetDataSet().GetDataElement( squenceTag ); + gdcm::DataSet &sequenceDataSet = sequence.GetValueAsSQ()->GetItem(1).GetNestedDataSet(); + gdcm::Tag rotationNumber_Tag(0x54,0x53); + const gdcm::DataElement &rotationNumber_Element = sequenceDataSet.GetDataElement( rotationNumber_Tag ); + gdcm::Attribute<0x54,0x53> rotationNumber_att; + rotationNumber_att.SetFromDataElement(rotationNumber_Element); + rotationNumber = rotationNumber_att.GetValue(); + +energyNumber = 3; +headNumber = 4; +rotationNumber = 6; + // Create a new mhd image with the dicom order slices + typename InputImageType::Pointer mhdCorrectOrder = InputImageType::New(); + mhdCorrectOrder->SetRegions(input->GetLargestPossibleRegion()); + mhdCorrectOrder->Allocate(); + unsigned int zAxis(0); //z value for the input mhd image + for (unsigned int energy = 0; energy < energyNumber; ++energy) { + for (unsigned int head = 0; head < headNumber; ++head) { + for (unsigned int rotation = 0; rotation < rotationNumber; ++rotation) { + std::cout << "Energy " << energy << " Head " << head << " Rotation " << rotation << std::endl; + + typename InputImageType::IndexType startIteratorIndexCorrectOrder; //pixel index of mhdCorrectOrder + startIteratorIndexCorrectOrder[0] = 0; + startIteratorIndexCorrectOrder[1] = 0; + startIteratorIndexCorrectOrder[2] = rotation + head*rotationNumber + energy*headNumber; + + typename InputImageType::IndexType startIteratorIndexOriginalOrder; //pixel index of input mhd + startIteratorIndexOriginalOrder[0] = 0; + startIteratorIndexOriginalOrder[1] = 0; + startIteratorIndexOriginalOrder[2] = head + energy*headNumber + rotation*energyNumber; + + typename InputImageType::SizeType regionSizeIterator; + regionSizeIterator[0] = input->GetLargestPossibleRegion().GetSize()[0]; + regionSizeIterator[1] = input->GetLargestPossibleRegion().GetSize()[1]; + regionSizeIterator[2] = 1; + + typename InputImageType::RegionType regionIteratorCorrectOrder; + regionIteratorCorrectOrder.SetSize(regionSizeIterator); + regionIteratorCorrectOrder.SetIndex(startIteratorIndexCorrectOrder); + + typename InputImageType::RegionType regionIteratorOriginalOrder; + regionIteratorOriginalOrder.SetSize(regionSizeIterator); + regionIteratorOriginalOrder.SetIndex(startIteratorIndexOriginalOrder); + + itk::ImageRegionIterator CorrectOrderIterator(mhdCorrectOrder,regionIteratorCorrectOrder); + itk::ImageRegionIterator OriginalOrderIterator(input,regionIteratorOriginalOrder); + while(!CorrectOrderIterator.IsAtEnd()) { + CorrectOrderIterator.Set(OriginalOrderIterator.Get()); + ++CorrectOrderIterator; + ++OriginalOrderIterator; + } + + ++zAxis; + } + } + } - -/* // update output dicom keys/tags - //std::string seriesUIDkey = "0020|000e"; - for (unsigned int i = 0; i < numberOfKeysGiven; i++) { - std::string entryId(m_ArgsInfo.key_arg[i] ); - std::string value( m_ArgsInfo.tag_arg[i] ); - - itk::EncapsulateMetaData(reader->GetMetaDataDictionary(), entryId, value ); - } -*/ + typename ReaderSeriesType::DictionaryRawPointer inputDict = (*(readerSeries->GetMetaDataDictionaryArray()))[0]; + typename ReaderSeriesType::DictionaryArrayType outputArray; + typename ReaderSeriesType::DictionaryRawPointer outputDict = new typename ReaderSeriesType::DictionaryType; + CopyDictionary (*inputDict, *outputDict); + + std::string entryId, value; + entryId = "0054|0011"; + value = NumberToString(energyNumber); + itk::EncapsulateMetaData(*outputDict, "0054|0011", value ); + entryId = "0054|0021"; + value = NumberToString(headNumber); + itk::EncapsulateMetaData(*outputDict, "0054|0021", value ); + +outputArray.push_back(outputDict); // Output directory and filenames //itksys::SystemTools::MakeDirectory( m_ArgsInfo.outputFilename_arg ); // create if it doesn't exist - typedef itk::ImageFileWriter SeriesWriterType; - typename SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New(); + typedef itk::ImageFileWriter WriterType; + typedef itk::ImageSeriesWriter WriterSerieType; + typename WriterType::Pointer writer = WriterType::New(); + typename WriterSerieType::Pointer writerSerie = WriterSerieType::New(); - seriesWriter->SetInput( input ); - seriesWriter->SetImageIO( gdcmIO ); + writer->SetInput( mhdCorrectOrder ); + writer->SetImageIO( gdcmIO ); - seriesWriter->SetFileName( filename_out ); - seriesWriter->SetMetaDataDictionary(reader->GetMetaDataDictionary()); + writer->SetFileName( filename_out ); + //writer->SetMetaDataDictionary(outputDict); + writerSerie->SetInput( mhdCorrectOrder ); + writerSerie->SetImageIO( gdcmIO ); + typename ReaderSeriesType::FileNamesContainer fileNamesOutput; + fileNamesOutput.push_back(filename_out); + writerSerie->SetFileNames( fileNamesOutput ); + writerSerie->SetMetaDataDictionaryArray(&outputArray); // Write try { if (m_ArgsInfo.verbose_flag) - std::cout << seriesWriter << std::endl; - seriesWriter->Update(); + std::cout << writer << std::endl; + //writer->Update(); + writerSerie->Update(); } catch( itk::ExceptionObject & excp ) { std::cerr << "Error: Exception thrown while writing the series!!" << std::endl; std::cerr << excp << std::endl; } +#else + std::cout << "Use GDCM2" << std::endl; +#endif +} +void CopyDictionary (itk::MetaDataDictionary &fromDict, itk::MetaDataDictionary &toDict) +{ + typedef itk::MetaDataDictionary DictionaryType; + + DictionaryType::ConstIterator itr = fromDict.Begin(); + DictionaryType::ConstIterator end = fromDict.End(); + typedef itk::MetaDataObject< std::string > MetaDataStringType; + + while( itr != end ) + { + itk::MetaDataObjectBase::Pointer entry = itr->second; + + MetaDataStringType::Pointer entryvalue = + dynamic_cast( entry.GetPointer() ) ; + if( entryvalue ) + { + std::string tagkey = itr->first; + std::string tagvalue = entryvalue->GetMetaDataObjectValue(); + itk::EncapsulateMetaData(toDict, tagkey, tagvalue); + } + ++itr; + } } +template + std::string NumberToString ( T Number ) + { + std::ostringstream ss; + ss << Number; + return ss.str(); + } }//end clitk #endif //#define clitkGateSimulation2DicomGenericFilter_txx -- 2.47.1