X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;f=tools%2FclitkPartitionEnergyWindowDicomGenericFilter.txx;fp=tools%2FclitkPartitionEnergyWindowDicomGenericFilter.txx;h=a32dc96d81a1076b8a0ab5cf42537216cf7fe41a;hb=3622e36be67020b2b9bd2c72c712a7a33b59b6e9;hp=0000000000000000000000000000000000000000;hpb=4916c13c75203b29bf65116afd754370f077e669;p=clitk.git diff --git a/tools/clitkPartitionEnergyWindowDicomGenericFilter.txx b/tools/clitkPartitionEnergyWindowDicomGenericFilter.txx new file mode 100644 index 0000000..a32dc96 --- /dev/null +++ b/tools/clitkPartitionEnergyWindowDicomGenericFilter.txx @@ -0,0 +1,506 @@ +/*========================================================================= + Program: vv http://www.creatis.insa-lyon.fr/rio/vv + + Authors belong to: + - University of LYON http://www.universite-lyon.fr/ + - Léon Bérard cancer center http://www.centreleonberard.fr + - CREATIS CNRS laboratory http://www.creatis.insa-lyon.fr + + This software is distributed WITHOUT ANY WARRANTY; without even + the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR + PURPOSE. See the copyright notices for more information. + + It is distributed under dual licence + + - BSD See included LICENSE.txt file + - CeCILL-B http://www.cecill.info/licences/Licence_CeCILL-B_V1-en.html +===========================================================================**/ +#ifndef clitkPartitionEnergyWindowDicomGenericFilter_txx +#define clitkPartitionEnergyWindowDicomGenericFilter_txx + +/* ================================================= + * @file clitkPartitionEnergyWindowDicomGenericFilter.txx + * @author + * @date + * + * @brief + * + ===================================================*/ + +#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 +{ + + +//----------------------------------------------------------- +// Constructor +//----------------------------------------------------------- +template +PartitionEnergyWindowDicomGenericFilter::PartitionEnergyWindowDicomGenericFilter() +{ + m_Verbose=false; + m_InputFileName=""; +} + + +//----------------------------------------------------------- +// Update +//----------------------------------------------------------- +template +void PartitionEnergyWindowDicomGenericFilter::Update() +{ + // Read the Dimension and PixelType + int Dimension; + std::string PixelType; + ReadImageDimensionAndPixelType(m_InputFileName, Dimension, PixelType); + + + // Call UpdateWithDim + if(Dimension==2) UpdateWithDim<2>(PixelType); + else if(Dimension==3) UpdateWithDim<3>(PixelType); + // else if (Dimension==4)UpdateWithDim<4>(PixelType); + else { + std::cout<<"Error, Only for 2 or 3 Dimensions!!!"< +template +void +PartitionEnergyWindowDicomGenericFilter::UpdateWithDim(std::string PixelType) +{ + if (m_Verbose) std::cout << "Image was detected to be "<(); + } + else if(PixelType == "unsigned_short"){ + if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_short..." << std::endl; + UpdateWithDimAndPixelType(); + } + + else if (PixelType == "unsigned_char") { + if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and unsigned_char..." << std::endl; + UpdateWithDimAndPixelType(); + } + + // else if (PixelType == "char"){ + // if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and signed_char..." << std::endl; + // UpdateWithDimAndPixelType(); + // } + else if (PixelType == "double") { + if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and double..." << std::endl; + UpdateWithDimAndPixelType(); + } + else { + if (m_Verbose) std::cout << "Launching filter in "<< Dimension <<"D and float..." << std::endl; + UpdateWithDimAndPixelType(); + } +} + +//------------------------------------------------------------------- +// Update with the number of dimensions and the pixeltype read from +// the dicom files. The MHD files may be resampled to match the +// dicom spacing (and number of slices). Rounding errors in resampling +// are handled by removing files when generating the output dicom +// series. +//------------------------------------------------------------------- +template +template +void +PartitionEnergyWindowDicomGenericFilter::UpdateWithDimAndPixelType() +{ + +#if GDCM_MAJOR_VERSION == 2 + // ImageTypes + typedef itk::Image InputImageType; + typedef itk::Image OutputImageType; + typedef itk::ImageSeriesReader< InputImageType > ReaderType; + typedef itk::ImageSeriesWriter< InputImageType, InputImageType > WriterType; + typedef itk::GDCMImageIO ImageIOType; + typedef typename InputImageType::RegionType RegionType; + typedef typename RegionType::SizeType SizeType; + + + // 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 ); + typename ReaderType::FileNamesContainer fileNamesInput; + fileNamesInput.push_back(m_ArgsInfo.input_arg); + reader->SetFileNames( fileNamesInput ); + try { + reader->Update(); + } catch (itk::ExceptionObject &excp) { + std::cerr << "Error: Exception thrown while reading the DICOM model file !!" << std::endl; + std::cerr << excp << std::endl; + } + typename InputImageType::Pointer input = reader->GetOutput(); + + typename ReaderType::DictionaryRawPointer inputDict = (*(reader->GetMetaDataDictionaryArray()))[0]; + 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)); + outputArray.push_back(outputDict); + + + //Read the number of slices and energy window + int nbSlice = input->GetLargestPossibleRegion().GetSize()[2]; + gdcm::Reader hreader; + hreader.SetFileName(m_ArgsInfo.input_arg); + hreader.Read(); + gdcm::DataSet& ds = hreader.GetFile().GetDataSet(); + gdcm::Attribute<0x54,0x11> series_number_att; + series_number_att.SetFromDataSet(ds); + int nbEnergyWindow = series_number_att.GetValue(); + + std::cout << nbSlice << " " << nbEnergyWindow << std::endl; + + int nbSlicePerEnergy = nbSlice / nbEnergyWindow; + if (nbSlicePerEnergy*nbEnergyWindow != nbSlice) + std::cerr << "Error: The number of slices is not correct !!" << std::endl; + + for (unsigned int energy = 0; energy < nbEnergyWindow; ++energy) { + + //Create the output + typename InputImageType::Pointer energyImage = InputImageType::New(); + energyImage->SetRegions(input->GetLargestPossibleRegion()); + typename InputImageType::IndexType startIndex; + startIndex[0] = 0; + startIndex[1] = 0; + startIndex[2] = 0;//energy*nbSlicePerEnergy; + typename InputImageType::SizeType regionSize; + regionSize[0] = input->GetLargestPossibleRegion().GetSize()[0]; + regionSize[1] = input->GetLargestPossibleRegion().GetSize()[1]; + regionSize[2] = nbSlicePerEnergy; + typename InputImageType::RegionType region; + region.SetSize(regionSize); + region.SetIndex(startIndex); + energyImage->SetRegions(region); + energyImage->Allocate(); + + //Create the iterator on the output + itk::ImageRegionIterator imageOutputIterator(energyImage,energyImage->GetLargestPossibleRegion()); + + //Create the iterator on the region of the input + typename InputImageType::IndexType startIndexIterator; + startIndexIterator[0] = 0; + startIndexIterator[1] = 0; + startIndexIterator[2] = energy*nbSlicePerEnergy; + typename InputImageType::RegionType regionIterator; + regionIterator.SetSize(regionSize); + regionIterator.SetIndex(startIndexIterator); + itk::ImageRegionIterator imageInputIterator(input,regionIterator); + + //Copy the requested region + while(!imageInputIterator.IsAtEnd()) + { + imageOutputIterator.Set(imageInputIterator.Get()); + + ++imageInputIterator; + ++imageOutputIterator; + } + + // Output directory and filenames + typename WriterType::Pointer writer = WriterType::New(); + writer->SetInput( energyImage ); + writer->SetImageIO( gdcmIO ); + typename ReaderType::FileNamesContainer fileNamesOutput; + fileNamesOutput.push_back(filename_out); + writer->SetFileNames( fileNamesOutput ); + writer->SetMetaDataDictionaryArray(&outputArray); + + // Write + try { + if (m_ArgsInfo.verbose_flag) + std::cout << writer << std::endl; + writer->Update(); + } catch( itk::ExceptionObject & excp ) { + 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(); + // 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 + +#endif //#define clitkPartitionEnergyWindowDicomGenericFilter_txx