Read, manipulate tags and write dicom series

I am reading dicom series and manipulate some tags from it. Finally write the new series. I can change any dicom tag without any problems. They are correcly changed. But the new dataset is visualized wrongly. The axial slice is not ordered correctly and many of the slices are missing. For sagittal and coronal views, I get only white lines and dots. I would be grateful for any suggestion, where the issue could be.

const unsigned int Dimension = 3;
const unsigned int OutputDimension = 2;

typedef signed short PixelType;
typedef itk::Image< PixelType, Dimension > InputImageType;
typedef itk::Image< PixelType, OutputDimension > OutputImageType;
typedef itk::GDCMSeriesFileNames InputNamesGeneratorType;
typedef itk::NumericSeriesFileNames OutputNamesGeneratorType;
typedef itk::ImageSeriesReader< InputImageType > ReaderType;
typedef itk::ImageSeriesWriter< InputImageType, OutputImageType >  SeriesWriterType;
typedef itk::GDCMImageIO ImageIOType;

ImageIOType::Pointer gdcmIO = ImageIOType::New();
InputNamesGeneratorType::Pointer inputNames = InputNamesGeneratorType::New();
inputNames->SetInputDirectory("C:/dicomDir/");
const ReaderType::FileNamesContainer& filenames = inputNames->GetInputFileNames();

// Read
ReaderType::Pointer seriesReader = ReaderType::New();
seriesReader->SetImageIO(gdcmIO);
seriesReader->SetFileNames(filenames);
seriesReader->Update();

//const InputImageType::SpacingType& inputSpacing = seriesReader->GetOutput()->GetSpacing();
const InputImageType::RegionType& inputRegion = seriesReader->GetOutput()->GetLargestPossibleRegion();
const InputImageType::SizeType& inputSize = inputRegion.GetSize();

InputImageType::SizeType   outputSize;
typedef InputImageType::SizeType::SizeValueType SizeValueType;
outputSize[0] = inputSize[0];
outputSize[1] = inputSize[1];
outputSize[2] = inputSize[2];

ReaderType::DictionaryRawPointer inputDict = (*(seriesReader->GetMetaDataDictionaryArray()))[0];
ReaderType::DictionaryArrayType outputArray;

gdcm::UIDGenerator sUid;
std::string seriesUid = sUid.Generate();
gdcm::UIDGenerator frameUid;
std::string frameOfReferenceUid = frameUid.Generate();

std::string studyUid;
std::string sopClassUid;
itk::ExposeMetaData<std::string>(*inputDict, "0020|000d", studyUid);
itk::ExposeMetaData<std::string>(*inputDict, "0008|0016", sopClassUid);
gdcmIO->KeepOriginalUIDOn();

for (unsigned int i = 0; i < outputSize[2]; i++)
{
  // Create a new dictionary for each slice
  ReaderType::DictionaryRawPointer dict = new ReaderType::DictionaryType;

  // Copy the dictionary from the first slice
  copyDictionary(*inputDict, *dict);

  // Set the UID's for the study, series, SOP and frame of reference
  itk::EncapsulateMetaData<std::string>(*dict, "0020|000d", studyUid); // Study Instance UID
  itk::EncapsulateMetaData<std::string>(*dict, "0020|000e", seriesUid); // Series Instance UID
  itk::EncapsulateMetaData<std::string>(*dict, "0020|0052", frameOfReferenceUid); // Frame of Reference UID

  gdcm::UIDGenerator sopUid;
  std::string sopInstanceUid = sopUid.Generate();
  itk::EncapsulateMetaData<std::string>(*dict, "0008|0018", sopInstanceUid); // SOP Instance UID
  itk::EncapsulateMetaData<std::string>(*dict, "0002|0003", sopInstanceUid); // Media Stored SOP Instance UID

    // Change patient name
    itk::EncapsulateMetaData<std::string>(*dict, "0010|0010", "Patient Name Anonymized"); // Patient Name

   // Change fields that are slice specific
    itksys_ios::ostringstream value;
    value.str("");
    value << i + 1;

    // Image Number
    itk::EncapsulateMetaData<std::string>(*dict, "0020|0013", value.str());
    outputArray.push_back(dict);
}

// Write
itksys::SystemTools::MakeDirectory("anonymizedPatientDataset");

OutputNamesGeneratorType::Pointer outputNames = OutputNamesGeneratorType::New();
std::string seriesFormat("anonymizedPatientDataset");
seriesFormat = seriesFormat + "/" + "000%d.dcm";
outputNames->SetSeriesFormat(seriesFormat.c_str());
outputNames->SetStartIndex(1);
outputNames->SetEndIndex(outputSize[2]);

SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New();
seriesWriter->SetInput(seriesReader->GetOutput());
seriesWriter->SetImageIO(gdcmIO);
seriesWriter->SetFileNames(outputNames->GetFileNames());
seriesWriter->SetMetaDataDictionaryArray(&outputArray);
seriesWriter->Update();

If you write seriesReader->GetOutput() as a 3D image, does it look correct? If not, your problem is in reading. Example: itk::WriteImage(seriesReader->GetOutput(), "test.nrrd");.

It might be useful taking a look at ResampleDICOMSeries example, if you haven’t already.

@csad6517 in 3D Slicer make sure you load the DICOM data set using the DICOM module. If you load using “Add data” then the slices may not be sorted.

Thank you for your suggestions. I don’t think that the issue neither in reading, resampling nor writing the dataset. I have integrated a resampler, the result does not change. Strangely, if I test the same code with a higher itk version e.g. 4.11. instead of my current version 3.20., the new series are written/visualized correctly.

const unsigned int Dimension = 3;
const unsigned int OutputDimension = 2;

typedef float PixelType;
typedef itk::Image< PixelType, Dimension > InputImageType;
typedef itk::Image< PixelType, OutputDimension > OutputImageType;
typedef itk::GDCMSeriesFileNames InputNamesGeneratorType;
typedef itk::NumericSeriesFileNames OutputNamesGeneratorType;
typedef itk::ImageSeriesReader< InputImageType > ReaderType;
typedef itk::ImageSeriesWriter< InputImageType, OutputImageType >  SeriesWriterType;
typedef itk::GDCMImageIO ImageIOType;
typedef itk::ResampleImageFilter<InputImageType, InputImageType> ResampleFilterType;

ImageIOType::Pointer gdcmIO = ImageIOType::New();
InputNamesGeneratorType::Pointer inputNames = InputNamesGeneratorType::New();
inputNames->SetInputDirectory("C:/dicomDir/");
const ReaderType::FileNamesContainer& filenames = inputNames->GetInputFileNames();

// Read
ReaderType::Pointer seriesReader = ReaderType::New();
seriesReader->SetImageIO(gdcmIO);
seriesReader->SetFileNames(filenames);
seriesReader->Update();

const InputImageType::RegionType &inputRegion = seriesReader->GetOutput()->GetLargestPossibleRegion();
const InputImageType::SizeType &inputSize = inputRegion.GetSize();

// Image size, origin, direction and spacing is the same as the original dataset
ResampleFilterType::Pointer resampler = ResampleFilterType::New();
resampler->SetInput(seriesReader->GetOutput());
resampler->SetOutputOrigin(seriesReader->GetOutput()->GetOrigin());
resampler->SetOutputSpacing(seriesReader->GetOutput()->GetSpacing());
resampler->SetOutputDirection(seriesReader->GetOutput()->GetDirection());
resampler->SetSize(inputSize);
resampler->Update();

ReaderType::DictionaryRawPointer inputDict = (*(seriesReader->GetMetaDataDictionaryArray()))[0];
ReaderType::DictionaryArrayType outputArray;

// Keep the new series and study UID's in the same study as the original, generate new series and frame of reference UID's.
gdcm::UIDGenerator serUid;
std::string seriesUid = serUid.Generate();
gdcm::UIDGenerator frameUid;
std::string frameOfReferenceUid = frameUid.Generate();

std::string studyUid;
std::string sopClassUid;
itk::ExposeMetaData<std::string>(*inputDict, "0020|000d", studyUid);
itk::ExposeMetaData<std::string>(*inputDict, "0008|0016", sopClassUid);
gdcmIo->KeepOriginalUIDOn();

for (unsigned int i = 0; i < inputSize[2]; i++)
{
  // Create a new dictionary for each slice
  ReaderType::DictionaryRawPointer dict = new ReaderType::DictionaryType;

  // Copy the dictionary from the first slice
  copyDictionary(*inputDict, *dict);

  // Set the UID's for the study, series, SOP and frame of reference
  itk::EncapsulateMetaData<std::string>(*dict, "0020|000d", studyUid); // Study Instance UID
  itk::EncapsulateMetaData<std::string>(*dict, "0020|000e", seriesUid); // Series Instance UID
  itk::EncapsulateMetaData<std::string>(*dict, "0020|0052", frameOfReferenceUid); // Frame of Reference UID

  gdcm::UIDGenerator sopUid;
  std::string sopInstanceUid = sopUid.Generate();
  itk::EncapsulateMetaData<std::string>(*dict, "0008|0018", sopInstanceUid); // SOP Instance UID
  itk::EncapsulateMetaData<std::string>(*dict, "0002|0003", sopInstanceUid); // Media Stored SOP Instance UID

  // Change patient name
  itk::EncapsulateMetaData<std::string>(*dict, "0010|0010", "Patient Name Anonymized"); // Patient Name

  // Change fields that are slice specific
  // Image Number
  std::ostringstream value;
  value.str("");
  value << i + 1;
  itk::EncapsulateMetaData<std::string>(*dict, "0020|0013", value.str());

  // Image Position Patient: This is calculated by computing the physical coordinate of the first pixel in each slice.
  InputImageType::PointType position;
  InputImageType::IndexType index;
  index[0] = 0;
  index[1] = 0;
  index[2] = i;
  resampler->GetOutput()->TransformIndexToPhysicalPoint(index, position);

  value.str("");
  value << position[0] << "\\" << position[1] << "\\" << position[2];
  itk::EncapsulateMetaData<std::string>(*dict, "0020|0032", value.str());

  // Slice Location: Store the z component of the Image Position Patient.
  value.str("");
  value << position[2];
  itk::EncapsulateMetaData<std::string>(*dict, "0020|1041", value.str());

  outputArray.push_back(dict);
}

// Write
itksys::SystemTools::MakeDirectory("anonymizedPatientDataset");

OutputNamesGeneratorType::Pointer outputNames = OutputNamesGeneratorType::New();
std::string seriesFormat("anonymizedPatientDataset");
seriesFormat = seriesFormat + "/" + "000%d.dcm";
outputNames->SetSeriesFormat(seriesFormat.c_str());
outputNames->SetStartIndex(1);
outputNames->SetEndIndex(inputSize[2]);

SeriesWriterType::Pointer seriesWriter = SeriesWriterType::New();
seriesWriter->SetInput(resampler->GetOutput());
seriesWriter->SetImageIO(gdcmIO);
seriesWriter->SetFileNames(outputNames->GetFileNames());
seriesWriter->SetMetaDataDictionaryArray(&outputArray);
seriesWriter->Update();

ITK 3.20 was released in 2010. You are missing on a decade of bug fixes.

1 Like

I realized that the 3th value of Image Position Patient Tag 0020|0032 is not ident with the Slice Location Tag 0020|1041 if I reslice the images, which is done exactly here:

	    // Calculated by computing the physical coordinate of the first pixel in each slice.
	    InputImageType::PointType position;
	    InputImageType::IndexType index;
	    index[0] = 0;
	    index[1] = 0;
	    index[2] = i;
	    resampler->GetOutput()->TransformIndexToPhysicalPoint(index, position);
        cout << "INDEX: " << index << "POSITION: " << position << endl;

	    value.str("");
	    value << position[0] << "//" << position[1] << "//" << position[2];
	    itk::EncapsulateMetaData<std::string>(*outputDict, "0020|0032", value.str()); // Image Position Patient

	    // Store the z component of the Image Position Patient.
	    value.str("");
	    value << position[2];
	    itk::EncapsulateMetaData<std::string>(*outputDict, "0020|1041", value.str()); // Slice Location

	    // Save the dictionary
	    outputArray.push_back(outputDict);

The Slice Location values are written correctly for each slice but the Image Position Patient values are written only for the last slice and it is constant (347.5 from index or rather slice number 157):

INDEX: [0, 0, 65]POSITION: [-109.781, -276.781, 283.1]
INDEX: [0, 0, 66]POSITION: [-109.781, -276.781, 283.8]
INDEX: [0, 0, 67]POSITION: [-109.781, -276.781, 284.5]
INDEX: [0, 0, 68]POSITION: [-109.781, -276.781, 285.2]
INDEX: [0, 0, 69]POSITION: [-109.781, -276.781, 285.9]
INDEX: [0, 0, 70]POSITION: [-109.781, -276.781, 286.6]
INDEX: [0, 0, 71]POSITION: [-109.781, -276.781, 287.3]
INDEX: [0, 0, 72]POSITION: [-109.781, -276.781, 288]
INDEX: [0, 0, 73]POSITION: [-109.781, -276.781, 288.7]
INDEX: [0, 0, 74]POSITION: [-109.781, -276.781, 289.4]
INDEX: [0, 0, 75]POSITION: [-109.781, -276.781, 290.1]
INDEX: [0, 0, 76]POSITION: [-109.781, -276.781, 290.8]
INDEX: [0, 0, 77]POSITION: [-109.781, -276.781, 291.5]
INDEX: [0, 0, 78]POSITION: [-109.781, -276.781, 292.2]
INDEX: [0, 0, 79]POSITION: [-109.781, -276.781, 292.9]
INDEX: [0, 0, 80]POSITION: [-109.781, -276.781, 293.6]
INDEX: [0, 0, 81]POSITION: [-109.781, -276.781, 294.3]
INDEX: [0, 0, 82]POSITION: [-109.781, -276.781, 295]
INDEX: [0, 0, 83]POSITION: [-109.781, -276.781, 295.7]
INDEX: [0, 0, 84]POSITION: [-109.781, -276.781, 296.4]
INDEX: [0, 0, 85]POSITION: [-109.781, -276.781, 297.1]
INDEX: [0, 0, 86]POSITION: [-109.781, -276.781, 297.8]
INDEX: [0, 0, 87]POSITION: [-109.781, -276.781, 298.5]
INDEX: [0, 0, 88]POSITION: [-109.781, -276.781, 299.2]
INDEX: [0, 0, 89]POSITION: [-109.781, -276.781, 299.9]
INDEX: [0, 0, 90]POSITION: [-109.781, -276.781, 300.6]
INDEX: [0, 0, 91]POSITION: [-109.781, -276.781, 301.3]
INDEX: [0, 0, 92]POSITION: [-109.781, -276.781, 302]
INDEX: [0, 0, 93]POSITION: [-109.781, -276.781, 302.7]
INDEX: [0, 0, 94]POSITION: [-109.781, -276.781, 303.4]
INDEX: [0, 0, 95]POSITION: [-109.781, -276.781, 304.1]
INDEX: [0, 0, 96]POSITION: [-109.781, -276.781, 304.8]
INDEX: [0, 0, 97]POSITION: [-109.781, -276.781, 305.5]
INDEX: [0, 0, 98]POSITION: [-109.781, -276.781, 306.2]
INDEX: [0, 0, 99]POSITION: [-109.781, -276.781, 306.9]
INDEX: [0, 0, 100]POSITION: [-109.781, -276.781, 307.6]
INDEX: [0, 0, 101]POSITION: [-109.781, -276.781, 308.3]
INDEX: [0, 0, 102]POSITION: [-109.781, -276.781, 309]
INDEX: [0, 0, 103]POSITION: [-109.781, -276.781, 309.7]
INDEX: [0, 0, 104]POSITION: [-109.781, -276.781, 310.4]
INDEX: [0, 0, 105]POSITION: [-109.781, -276.781, 311.1]
INDEX: [0, 0, 106]POSITION: [-109.781, -276.781, 311.8]
INDEX: [0, 0, 107]POSITION: [-109.781, -276.781, 312.5]
INDEX: [0, 0, 108]POSITION: [-109.781, -276.781, 313.2]
INDEX: [0, 0, 109]POSITION: [-109.781, -276.781, 313.9]
INDEX: [0, 0, 110]POSITION: [-109.781, -276.781, 314.6]
INDEX: [0, 0, 111]POSITION: [-109.781, -276.781, 315.3]
INDEX: [0, 0, 112]POSITION: [-109.781, -276.781, 316]
INDEX: [0, 0, 113]POSITION: [-109.781, -276.781, 316.7]
INDEX: [0, 0, 114]POSITION: [-109.781, -276.781, 317.4]
INDEX: [0, 0, 115]POSITION: [-109.781, -276.781, 318.1]
INDEX: [0, 0, 116]POSITION: [-109.781, -276.781, 318.8]
INDEX: [0, 0, 117]POSITION: [-109.781, -276.781, 319.5]
INDEX: [0, 0, 118]POSITION: [-109.781, -276.781, 320.2]
INDEX: [0, 0, 119]POSITION: [-109.781, -276.781, 320.9]
INDEX: [0, 0, 120]POSITION: [-109.781, -276.781, 321.6]
INDEX: [0, 0, 121]POSITION: [-109.781, -276.781, 322.3]
INDEX: [0, 0, 122]POSITION: [-109.781, -276.781, 323]
INDEX: [0, 0, 123]POSITION: [-109.781, -276.781, 323.7]
INDEX: [0, 0, 124]POSITION: [-109.781, -276.781, 324.4]
INDEX: [0, 0, 125]POSITION: [-109.781, -276.781, 325.1]
INDEX: [0, 0, 126]POSITION: [-109.781, -276.781, 325.8]
INDEX: [0, 0, 127]POSITION: [-109.781, -276.781, 326.5]
INDEX: [0, 0, 128]POSITION: [-109.781, -276.781, 327.2]
INDEX: [0, 0, 129]POSITION: [-109.781, -276.781, 327.9]
INDEX: [0, 0, 130]POSITION: [-109.781, -276.781, 328.6]
INDEX: [0, 0, 131]POSITION: [-109.781, -276.781, 329.3]
INDEX: [0, 0, 132]POSITION: [-109.781, -276.781, 330]
INDEX: [0, 0, 133]POSITION: [-109.781, -276.781, 330.7]
INDEX: [0, 0, 134]POSITION: [-109.781, -276.781, 331.4]
INDEX: [0, 0, 135]POSITION: [-109.781, -276.781, 332.1]
INDEX: [0, 0, 136]POSITION: [-109.781, -276.781, 332.8]
INDEX: [0, 0, 137]POSITION: [-109.781, -276.781, 333.5]
INDEX: [0, 0, 138]POSITION: [-109.781, -276.781, 334.2]
INDEX: [0, 0, 139]POSITION: [-109.781, -276.781, 334.9]
INDEX: [0, 0, 140]POSITION: [-109.781, -276.781, 335.6]
INDEX: [0, 0, 141]POSITION: [-109.781, -276.781, 336.3]
INDEX: [0, 0, 142]POSITION: [-109.781, -276.781, 337]
INDEX: [0, 0, 143]POSITION: [-109.781, -276.781, 337.7]
INDEX: [0, 0, 144]POSITION: [-109.781, -276.781, 338.4]
INDEX: [0, 0, 145]POSITION: [-109.781, -276.781, 339.1]
INDEX: [0, 0, 146]POSITION: [-109.781, -276.781, 339.8]
INDEX: [0, 0, 147]POSITION: [-109.781, -276.781, 340.5]
INDEX: [0, 0, 148]POSITION: [-109.781, -276.781, 341.2]
INDEX: [0, 0, 149]POSITION: [-109.781, -276.781, 341.9]
INDEX: [0, 0, 150]POSITION: [-109.781, -276.781, 342.6]
INDEX: [0, 0, 151]POSITION: [-109.781, -276.781, 343.3]
INDEX: [0, 0, 152]POSITION: [-109.781, -276.781, 344]
INDEX: [0, 0, 153]POSITION: [-109.781, -276.781, 344.7]
INDEX: [0, 0, 154]POSITION: [-109.781, -276.781, 345.4]
INDEX: [0, 0, 155]POSITION: [-109.781, -276.781, 346.1]
INDEX: [0, 0, 156]POSITION: [-109.781, -276.781, 346.8]
INDEX: [0, 0, 157]POSITION: [-109.781, -276.781, 347.5]

This was a known bug in 3.20.0 from 2011: https://itk.org/pipermail/insight-users/2011-August/042224.html . Unfortunately I can not follow the bug tracking, since the website is down