X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=tools%2FclitkPartitionEnergyWindowDicomGenericFilter.txx;h=8e8fa33f37f23359d785f9d6531b8e4aab316382;hb=5578995d9a82792833333eeb3dd5c8ecac967293;hp=a32dc96d81a1076b8a0ab5cf42537216cf7fe41a;hpb=3622e36be67020b2b9bd2c72c712a7a33b59b6e9;p=clitk.git diff --git a/tools/clitkPartitionEnergyWindowDicomGenericFilter.txx b/tools/clitkPartitionEnergyWindowDicomGenericFilter.txx index a32dc96..8e8fa33 100644 --- a/tools/clitkPartitionEnergyWindowDicomGenericFilter.txx +++ b/tools/clitkPartitionEnergyWindowDicomGenericFilter.txx @@ -36,6 +36,8 @@ #include #include #include +#include +#include #else #include "gdcmFile.h" #include "gdcmUtil.h" @@ -134,7 +136,7 @@ void PartitionEnergyWindowDicomGenericFilter::UpdateWithDimAndPixelType() { -#if GDCM_MAJOR_VERSION == 2 +#if GDCM_MAJOR_VERSION >= 2 // ImageTypes typedef itk::Image InputImageType; typedef itk::Image OutputImageType; @@ -148,7 +150,6 @@ PartitionEnergyWindowDicomGenericFilter::UpdateWithDimAndPixelTy // Read Dicom model file typename ReaderType::Pointer reader = ReaderType::New(); ImageIOType::Pointer gdcmIO = ImageIOType::New(); - std::string filename_out = m_ArgsInfo.output_arg; gdcmIO->LoadPrivateTagsOn(); gdcmIO->KeepOriginalUIDOn(); reader->SetImageIO( gdcmIO ); @@ -167,22 +168,33 @@ PartitionEnergyWindowDicomGenericFilter::UpdateWithDimAndPixelTy typename ReaderType::DictionaryArrayType outputArray; typename ReaderType::DictionaryRawPointer outputDict = new typename ReaderType::DictionaryType; CopyDictionary (*inputDict, *outputDict); - //itk::EncapsulateMetaData(*outputDict, "0054|0011", NumberToString(energyNumber)); - //itk::EncapsulateMetaData(*outputDict, "0054|0021", NumberToString(headNumber)); + itk::EncapsulateMetaData(*outputDict, "0054|0011", NumberToString(1)); outputArray.push_back(outputDict); - //Read the number of slices and energy window + //Read the number of slices and energy window, rotation int nbSlice = input->GetLargestPossibleRegion().GetSize()[2]; gdcm::Reader hreader; hreader.SetFileName(m_ArgsInfo.input_arg); hreader.Read(); - gdcm::DataSet& ds = hreader.GetFile().GetDataSet(); + gdcm::DataSet& dsInput = hreader.GetFile().GetDataSet(); gdcm::Attribute<0x54,0x11> series_number_att; - series_number_att.SetFromDataSet(ds); + series_number_att.SetFromDataSet(dsInput); int nbEnergyWindow = series_number_att.GetValue(); - std::cout << nbSlice << " " << nbEnergyWindow << std::endl; + int nbRotation; + if (dsInput.FindDataElement(gdcm::Tag(0x54, 0x52))) { + const gdcm::DataElement &seqRotation = dsInput.GetDataElement(gdcm::Tag(0x54, 0x52)); + gdcm::SmartPointer sqiRotation = seqRotation.GetValueAsSQ(); + gdcm::Item &itemRotation = sqiRotation->GetItem(1); + const gdcm::DataElement &deRotation = itemRotation.GetDataElement(gdcm::Tag(0x54, 0x53)); + unsigned short SamplesPerPixelRotation = *((unsigned short *) deRotation.GetByteValue()->GetPointer()); + std::stringstream s_SamplesPerPixelRotation; + s_SamplesPerPixelRotation << SamplesPerPixelRotation; + std::string p_StrRotation = s_SamplesPerPixelRotation.str(); + nbRotation = (int)atof(p_StrRotation.c_str()); + } else + nbRotation = 1; int nbSlicePerEnergy = nbSlice / nbEnergyWindow; if (nbSlicePerEnergy*nbEnergyWindow != nbSlice) @@ -234,6 +246,8 @@ PartitionEnergyWindowDicomGenericFilter::UpdateWithDimAndPixelTy writer->SetInput( energyImage ); writer->SetImageIO( gdcmIO ); typename ReaderType::FileNamesContainer fileNamesOutput; + std::string extension = NumberToString(energy) + ".dcm"; + std::string filename_out = m_ArgsInfo.output_arg + extension; fileNamesOutput.push_back(filename_out); writer->SetFileNames( fileNamesOutput ); writer->SetMetaDataDictionaryArray(&outputArray); @@ -247,217 +261,106 @@ PartitionEnergyWindowDicomGenericFilter::UpdateWithDimAndPixelTy std::cerr << "Error: Exception thrown while writing the series!!" << std::endl; std::cerr << excp << std::endl; } - } - -/* - - //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 - // 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); - - // Output directory and filenames - 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); - - - // Write - try { - if (m_ArgsInfo.verbose_flag) - 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(); + //Write sequence dicom tag with gdcm + gdcm::Reader reader; + reader.SetFileName( filename_out.c_str() ); + reader.Read(); + gdcm::File &fileOutput = reader.GetFile(); + gdcm::DataSet &dsOutput = fileOutput.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 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])); + gdcm::DataElement rotationDE( gdcm::Tag(0x54, 0x53) ); + rotationDE.SetVR( gdcm::VR::US ); + char essai = (char)nbRotation; + char *p_essai = &essai; + rotationDE.SetByteValue(p_essai, 1); // Create an item - gdcm::Item energyThresholdIt; - energyThresholdIt.SetVLToUndefined(); - gdcm::DataSet &energyThresholdDS = energyThresholdIt.GetNestedDataSet(); - energyThresholdDS.Insert(energyThresholdDE); - energyThresholdDS.Insert(energyUpholdDE); - energyThresholdSq->AddItem(energyThresholdIt); + gdcm::Item rotationIt; + rotationIt.SetVLToUndefined(); + gdcm::DataSet &rotationDS = rotationIt.GetNestedDataSet(); + rotationDS.Insert(rotationDE); + rotationSq->AddItem(rotationIt); // 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); + 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(); + + //Read the caracteristics of this energy window : the name and the thresholds up and down + const gdcm::DataElement &seqWindow = dsInput.GetDataElement(gdcm::Tag(0x54, 0x12)); + gdcm::SmartPointer sqiWindow = seqWindow.GetValueAsSQ(); + gdcm::Item &itemWindow = sqiWindow->GetItem(energy+1); + const gdcm::DataElement &deNameWindow = itemWindow.GetDataElement(gdcm::Tag(0x54, 0x18)); + std::string nameWindow( deNameWindow.GetByteValue()->GetPointer(), deNameWindow.GetByteValue()->GetLength() ); + + const gdcm::DataElement &deThresholds = itemWindow.GetDataElement(gdcm::Tag(0x54, 0x13)); + gdcm::SmartPointer sqiThresholds = deThresholds.GetValueAsSQ(); + for (unsigned int nbWindow=1; nbWindow <= sqiThresholds->GetNumberOfItems() ; ++nbWindow) { + gdcm::Item &itemThresholds = sqiThresholds->GetItem(nbWindow); + const gdcm::DataElement &deThresholdDown = itemThresholds.GetDataElement(gdcm::Tag(0x54, 0x14)); + std::string p_StrThresholdDown( deThresholdDown.GetByteValue()->GetPointer(), deThresholdDown.GetByteValue()->GetLength() ); + double thresholdDown = (double)atof(p_StrThresholdDown.c_str()); + const gdcm::DataElement &deThresholdUp = itemThresholds.GetDataElement(gdcm::Tag(0x54, 0x15)); + std::string p_StrThresholdUp( deThresholdUp.GetByteValue()->GetPointer(), deThresholdUp.GetByteValue()->GetLength() ); + double thresholdUp = (double)atof(p_StrThresholdUp.c_str()); + + //Now write it ! + gdcm::SmartPointer 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(NumberToString(thresholdDown).c_str(), (uint32_t)strlen(NumberToString(thresholdDown).c_str())); + energyUpholdDE.SetByteValue(NumberToString(thresholdUp).c_str(), (uint32_t)strlen(NumberToString(thresholdUp).c_str())); + // 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(nameWindow.c_str(), (uint32_t)strlen(nameWindow.c_str())); + // 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(); + dsOutput.Insert(energyDEParent); + dsOutput.Insert(rotationDEParent); + gdcm::Writer w; + w.SetFile( fileOutput ); + w.SetFileName( filename_out.c_str() ); + w.Write(); } - // 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