X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=tools%2FclitkGateSimulation2DicomGenericFilter.txx;h=9ad2857f90bd8bb263463b47b932f33cabaaed93;hb=85d87a4260fedf7ce53875ebd8654787a2dd941c;hp=f49bebcee94966dfe57ffb3588eee40882e77818;hpb=b155b373f6f32609a18c73e6756566947dd2f892;p=clitk.git diff --git a/tools/clitkGateSimulation2DicomGenericFilter.txx b/tools/clitkGateSimulation2DicomGenericFilter.txx index f49bebc..9ad2857 100644 --- a/tools/clitkGateSimulation2DicomGenericFilter.txx +++ b/tools/clitkGateSimulation2DicomGenericFilter.txx @@ -27,15 +27,24 @@ * ===================================================*/ +#include // clitk #include "clitkResampleImageWithOptionsFilter.h" #if GDCM_MAJOR_VERSION >= 2 #include "gdcmUIDGenerator.h" +#include +#include +#include +#include #else #include "gdcmFile.h" #include "gdcmUtil.h" #endif +#include "itkImageRegionIterator.h" +#include "itkMetaImageIO.h" +#include "itkMetaDataDictionary.h" + namespace clitk { @@ -125,39 +134,33 @@ 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; + + // Read Dicom model file + typename ReaderSeriesType::Pointer readerSeries = ReaderSeriesType::New(); 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 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,82 +168,258 @@ 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()) { - - // resampling is carried out on the fly if resolution or size between - // the input mhd and input dicom series is different - - // Filter - typedef clitk::ResampleImageWithOptionsFilter ResampleImageFilterType; - typename ResampleImageFilterType::Pointer filter = ResampleImageFilterType::New(); - filter->SetInput(input); - filter->SetVerboseOptions(m_Verbose); - filter->SetGaussianFilteringEnabled(false); - filter->SetDefaultPixelValue(0); - - const SizeType& input_size = input->GetLargestPossibleRegion().GetSize(); - SizeType output_size; - for (unsigned int i = 0; i < Dimension; i++) - output_size[i] = input_size[i]; - filter->SetOutputSize(output_size); - if (m_Verbose) { - std::cout << "Warning: The image size differs between the MHD file and the input dicom series. Performing resampling with default options using mhd size as reference (for advanced options, use clitkResampleImage)." << std::endl; - std::cout << "MHD -> " << input->GetLargestPossibleRegion().GetSize() << std::endl; - std::cout << "dicom -> " << reader->GetOutput()->GetLargestPossibleRegion().GetSize() << std::endl; - } - try { - filter->Update(); - input = filter->GetOutput(); - } catch( itk::ExceptionObject & excp ) { - std::cerr << "Error: Exception thrown while resampling!!" << std::endl; - 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) + //Read the metadata informations in the mhd + MetaObject tObj(3); + tObj.AddUserField("NumberOfEnergyWindows", MET_STRING); + tObj.AddUserField("NbOfHeads", MET_STRING); + tObj.AddUserField("NbOfProjections", MET_STRING); + tObj.Read(m_InputFileName.c_str()); + char *p_energyNumber, *p_headNumber, *p_rotationNumber; + int energyNumber, headNumber, rotationNumber; + p_energyNumber = static_cast(tObj.GetUserField("NumberOfEnergyWindows")); + p_headNumber = static_cast(tObj.GetUserField("NbOfHeads")); + p_rotationNumber = static_cast(tObj.GetUserField("NbOfProjections")); + energyNumber = atoi(p_energyNumber); + headNumber = atoi(p_headNumber); + rotationNumber = atoi(p_rotationNumber); + for (unsigned int i=0; i p_EnergyWindowName(energyNumber), p_EnergyWindowThreshold(energyNumber), p_EnergyWindowUphold(energyNumber); + std::vector energyWindowThreshold(energyNumber), energyWindowUphold(energyNumber); + for (unsigned int i=0; i(tObj.GetUserField(tagEnergyName.c_str())); + p_EnergyWindowThreshold[i] = static_cast(tObj.GetUserField(tagEnergyThreshold.c_str())); + p_EnergyWindowUphold[i] = static_cast(tObj.GetUserField(tagEnergyUphold.c_str())); + energyWindowThreshold[i] = atoi(p_EnergyWindowThreshold[i]); + energyWindowUphold[i] = atoi(p_EnergyWindowUphold[i]); + } +//TODO if the input size/spacing and dicom model ones are different + + // Create a new mhd image with the correct 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] ); + // string for distinguishing items inside sequence: + const std::string ITEM_ENCAPSULATE_STRING("DICOM_ITEM_ENCAPSULATE"); + std::string tempString = ITEM_ENCAPSULATE_STRING + "01"; + typename ReaderSeriesType::DictionaryRawPointer inputDict = (*(readerSeries->GetMetaDataDictionaryArray()))[0]; + typename ReaderSeriesType::DictionaryArrayType outputArray; + typename ReaderSeriesType::DictionaryRawPointer outputDict = new typename ReaderSeriesType::DictionaryType; + CopyDictionary (*inputDict, *outputDict); + itk::EncapsulateMetaData(*outputDict, "0054|0011", NumberToString(energyNumber)); + itk::EncapsulateMetaData(*outputDict, "0054|0021", NumberToString(headNumber)); + outputArray.push_back(outputDict); + - itk::EncapsulateMetaData(reader->GetMetaDataDictionary(), entryId, value ); - } -*/ // 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::ImageSeriesWriter WriterSerieType; + typename WriterSerieType::Pointer writerSerie = WriterSerieType::New(); + writerSerie->SetInput( mhdCorrectOrder ); + writerSerie->SetImageIO( gdcmIO ); + typename ReaderSeriesType::FileNamesContainer fileNamesOutput; + fileNamesOutput.push_back(filename_out); + writerSerie->SetFileNames( fileNamesOutput ); + writerSerie->SetMetaDataDictionaryArray(&outputArray); - seriesWriter->SetInput( input ); - seriesWriter->SetImageIO( gdcmIO ); - - seriesWriter->SetFileName( filename_out ); - seriesWriter->SetMetaDataDictionary(reader->GetMetaDataDictionary()); // Write try { if (m_ArgsInfo.verbose_flag) - std::cout << seriesWriter << std::endl; - seriesWriter->Update(); + std::cout << writerSerie << std::endl; + writerSerie->Update(); } catch( itk::ExceptionObject & excp ) { std::cerr << "Error: Exception thrown while writing the series!!" << std::endl; std::cerr << excp << std::endl; } + + //Write sequence dicom tag with gdcm + gdcm::Reader reader; + reader.SetFileName( fileNamesOutput[0].c_str() ); + reader.Read(); + gdcm::File &file = reader.GetFile(); + gdcm::DataSet &ds2 = file.GetDataSet(); + const unsigned int ptr_len = 42; + char *ptr = new char[ptr_len]; + memset(ptr,0,ptr_len); + + //Write rotation tag + // Create a Sequence + gdcm::SmartPointer rotationSq = new gdcm::SequenceOfItems(); + rotationSq->SetLengthToUndefined(); + // Create a dataelement + gdcm::DataElement rotationDE( gdcm::Tag(0x54, 0x53) ); + rotationDE.SetVR( gdcm::VR::US ); + char essai = (char)rotationNumber; + char *p_essai = &essai; + rotationDE.SetByteValue(p_essai, 1); + // Create an item + gdcm::Item rotationIt; + rotationIt.SetVLToUndefined(); + gdcm::DataSet &rotationDS = rotationIt.GetNestedDataSet(); + rotationDS.Insert(rotationDE); + rotationSq->AddItem(rotationIt); + // Insert sequence into data set + gdcm::DataElement rotationDEParent( gdcm::Tag(0x54, 0x52) ); + rotationDEParent.SetVR(gdcm::VR::SQ); + rotationDEParent.SetValue(*rotationSq); + rotationDEParent.SetVLToUndefined(); + + //Write energy + gdcm::DataElement energyDEParent( gdcm::Tag(0x54, 0x12) ); + energyDEParent.SetVR(gdcm::VR::SQ); + // Create a Sequence + gdcm::SmartPointer energySq = new gdcm::SequenceOfItems(); + energySq->SetLengthToUndefined(); + for (unsigned int i=0; i energyThresholdSq = new gdcm::SequenceOfItems(); + energyThresholdSq->SetLengthToUndefined(); + // Create a dataelement + gdcm::DataElement energyThresholdDE( gdcm::Tag(0x54, 0x14) ); + gdcm::DataElement energyUpholdDE( gdcm::Tag(0x54, 0x15) ); + energyThresholdDE.SetVR( gdcm::VR::DS ); + energyUpholdDE.SetVR( gdcm::VR::DS ); + energyThresholdDE.SetByteValue(p_EnergyWindowThreshold[i], (uint32_t)strlen(p_EnergyWindowThreshold[i])); + energyUpholdDE.SetByteValue(p_EnergyWindowUphold[i], (uint32_t)strlen(p_EnergyWindowUphold[i])); + // Create an item + gdcm::Item energyThresholdIt; + energyThresholdIt.SetVLToUndefined(); + gdcm::DataSet &energyThresholdDS = energyThresholdIt.GetNestedDataSet(); + energyThresholdDS.Insert(energyThresholdDE); + energyThresholdDS.Insert(energyUpholdDE); + energyThresholdSq->AddItem(energyThresholdIt); + // Insert sequence into data set + gdcm::DataElement energyThresholdDEParent( gdcm::Tag(0x54, 0x13) ); + energyThresholdDEParent.SetVR(gdcm::VR::SQ); + energyThresholdDEParent.SetValue(*energyThresholdSq); + energyThresholdDEParent.SetVLToUndefined(); + // Create a dataelement + gdcm::DataElement energyNameDE( gdcm::Tag(0x54, 0x18) ); + energyNameDE.SetVR( gdcm::VR::SH ); + energyNameDE.SetByteValue(p_EnergyWindowName[i], (uint32_t)strlen(p_EnergyWindowName[i])); + // Create an item + gdcm::Item energyIt; + energyIt.SetVLToUndefined(); + gdcm::DataSet &energyDS = energyIt.GetNestedDataSet(); + energyDS.Insert(energyNameDE); + energyDS.Insert(energyThresholdDEParent); + energySq->AddItem(energyIt); + } + // Insert sequence into data set + energyDEParent.SetValue(*energySq); + energyDEParent.SetVLToUndefined(); + ds2.Insert(energyDEParent); + ds2.Insert(rotationDEParent); + + gdcm::Writer w; + w.SetFile( file ); + w.SetFileName( fileNamesOutput[0].c_str() ); + w.Write(); + +#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