X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=tools%2FclitkGateSimulation2DicomGenericFilter.txx;h=9ad2857f90bd8bb263463b47b932f33cabaaed93;hb=342abd5c807dab47a30d8369aca8f7e5c9242b2b;hp=3b54bdca143063be7d1837b1cab90fdb41ce342d;hpb=15ebed734d7ef4727dff24b21d00d929b9162f90;p=clitk.git diff --git a/tools/clitkGateSimulation2DicomGenericFilter.txx b/tools/clitkGateSimulation2DicomGenericFilter.txx index 3b54bdc..9ad2857 100644 --- a/tools/clitkGateSimulation2DicomGenericFilter.txx +++ b/tools/clitkGateSimulation2DicomGenericFilter.txx @@ -42,6 +42,8 @@ #endif #include "itkImageRegionIterator.h" +#include "itkMetaImageIO.h" +#include "itkMetaDataDictionary.h" namespace clitk @@ -132,7 +134,7 @@ void GateSimulation2DicomGenericFilter::UpdateWithDimAndPixelType() { -#if GDCM_MAJOR_VERSION == 2 +#if GDCM_MAJOR_VERSION >= 2 // ImageTypes typedef itk::Image InputImageType; typedef itk::Image OutputImageType; @@ -142,20 +144,16 @@ GateSimulation2DicomGenericFilter::UpdateWithDimAndPixelType() // Read Dicom model file - typename ReaderType::Pointer reader = ReaderType::New(); typename ReaderSeriesType::Pointer readerSeries = ReaderSeriesType::New(); ImageIOType::Pointer gdcmIO = ImageIOType::New(); std::string filename_out = m_ArgsInfo.outputFilename_arg; gdcmIO->LoadPrivateTagsOn(); gdcmIO->KeepOriginalUIDOn(); - 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 model file !!" << std::endl; @@ -173,69 +171,49 @@ GateSimulation2DicomGenericFilter::UpdateWithDimAndPixelType() typename InputImageType::Pointer input = volumeReader->GetOutput(); + //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 -/* 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; - } - } -*/ - - - //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 + // Create a new mhd image with the correct dicom order slices typename InputImageType::Pointer mhdCorrectOrder = InputImageType::New(); mhdCorrectOrder->SetRegions(input->GetLargestPossibleRegion()); mhdCorrectOrder->Allocate(); @@ -289,34 +267,15 @@ rotationNumber = 6; typename ReaderSeriesType::DictionaryRawPointer inputDict = (*(readerSeries->GetMetaDataDictionaryArray()))[0]; typename ReaderSeriesType::DictionaryArrayType outputArray; typename ReaderSeriesType::DictionaryRawPointer outputDict = new typename ReaderSeriesType::DictionaryType; - typename ReaderSeriesType::DictionaryRawPointer sequenceDict = new typename ReaderSeriesType::DictionaryType; - // Dictionary for encapsulating the sequences. - itk::GDCMImageIO::Pointer gdcmIOStructureSetRoi = itk::GDCMImageIO::New(); - typename ReaderSeriesType::DictionaryType &structureSetRoiSeq = gdcmIOStructureSetRoi->GetMetaDataDictionary(); CopyDictionary (*inputDict, *outputDict); - //CopyDictionary (*inputDict, *sequenceDict); - itk::EncapsulateMetaData(*outputDict, "0054|0011", NumberToString(energyNumber)); itk::EncapsulateMetaData(*outputDict, "0054|0021", NumberToString(headNumber)); + outputArray.push_back(outputDict); - itk::EncapsulateMetaData(*sequenceDict, "0054|0053", NumberToString(rotationNumber)); - //itk::EncapsulateMetaData(*outputDict, "0054|0053", NumberToString(rotationNumber)); - itk::EncapsulateMetaData(structureSetRoiSeq, tempString, *sequenceDict); - itk::EncapsulateMetaData(*outputDict, "0054|0052",structureSetRoiSeq); -outputArray.push_back(outputDict); // Output directory and filenames - //itksys::SystemTools::MakeDirectory( m_ArgsInfo.outputFilename_arg ); // create if it doesn't exist - typedef itk::ImageFileWriter WriterType; typedef itk::ImageSeriesWriter WriterSerieType; - typename WriterType::Pointer writer = WriterType::New(); typename WriterSerieType::Pointer writerSerie = WriterSerieType::New(); - - writer->SetInput( mhdCorrectOrder ); - writer->SetImageIO( gdcmIO ); - - writer->SetFileName( filename_out ); - //writer->SetMetaDataDictionary(outputDict); writerSerie->SetInput( mhdCorrectOrder ); writerSerie->SetImageIO( gdcmIO ); typename ReaderSeriesType::FileNamesContainer fileNamesOutput; @@ -324,65 +283,98 @@ outputArray.push_back(outputDict); writerSerie->SetFileNames( fileNamesOutput ); writerSerie->SetMetaDataDictionaryArray(&outputArray); + // Write try { if (m_ArgsInfo.verbose_flag) - std::cout << writer << std::endl; - //writer->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; } - gdcm::Reader reader2; - reader2.SetFileName( fileNamesOutput[0].c_str() ); - reader2.Read(); - gdcm::File &file = reader2.GetFile(); + //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 nitems = 1000; - const unsigned int ptr_len = 42; /*94967296 / nitems; */ - //assert( ptr_len == 42949672 ); + 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 sq = new gdcm::SequenceOfItems(); - sq->SetLengthToUndefined(); - - //const char owner_str[] = "GDCM CONFORMANCE TESTS"; - //gdcm::DataElement owner( gdcm::Tag(0x54, 0x52) ); - //owner.SetByteValue(owner_str, (uint32_t)strlen(owner_str)); - //owner.SetVR( gdcm::VR::LO ); - + gdcm::SmartPointer rotationSq = new gdcm::SequenceOfItems(); + rotationSq->SetLengthToUndefined(); // Create a dataelement - gdcm::DataElement de( gdcm::Tag(0x54, 0x53) ); - de.SetByteValue(NumberToString(rotationNumber).c_str(), 1); - de.SetVR( gdcm::VR::OB ); - + 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 it; - it.SetVLToUndefined(); - gdcm::DataSet &nds = it.GetNestedDataSet(); - //nds.Insert(owner); - nds.Insert(de); - - sq->AddItem(it); - + gdcm::Item rotationIt; + rotationIt.SetVLToUndefined(); + gdcm::DataSet &rotationDS = rotationIt.GetNestedDataSet(); + rotationDS.Insert(rotationDE); + rotationSq->AddItem(rotationIt); // Insert sequence into data set - gdcm::DataElement des( gdcm::Tag(0x54, 0x52) ); - des.SetVR(gdcm::VR::SQ); - des.SetValue(*sq); - des.SetVLToUndefined(); - - //ds2.Insert(owner); - ds2.Insert(des); + 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.SetCheckFileMetaInformation( true ); w.SetFileName( fileNamesOutput[0].c_str() ); w.Write(); @@ -390,7 +382,10 @@ outputArray.push_back(outputDict); std::cout << "Use GDCM2" << std::endl; #endif } +//--------------------------------------------------------------------------- + +//--------------------------------------------------------------------------- void CopyDictionary (itk::MetaDataDictionary &fromDict, itk::MetaDataDictionary &toDict) { typedef itk::MetaDataDictionary DictionaryType; @@ -414,14 +409,18 @@ void CopyDictionary (itk::MetaDataDictionary &fromDict, itk::MetaDataDictionary ++itr; } } +//--------------------------------------------------------------------------- + + +//--------------------------------------------------------------------------- +template std::string NumberToString ( T Number ) +{ + std::ostringstream ss; + ss << Number; + return ss.str(); +} +//--------------------------------------------------------------------------- -template - std::string NumberToString ( T Number ) - { - std::ostringstream ss; - ss << Number; - return ss.str(); - } }//end clitk #endif //#define clitkGateSimulation2DicomGenericFilter_txx