Orientation in Diffusion Tensor Imaging file (DTI) .nii

Hello Everyone!

I have two .nii files, one containing a tumor segmentation, and one contains a DTI dataset. These two datasets have the same orientations/offsets in their metadata. I run numerical simulations using the DTI tensor data, and want to write back to a .nii file format so that my medical partners can make use of the results.

I have a problem with the orientations, which are stored in the metadata. I can read the files just fine, but when I write my results back to file, the orientations and offsets are all over the place.

I tried to just copy the metadata from the segmentation file, in the hope of retaining the orientation information, but this does not seem to work.

Reading segmentation and metadata:

	void read_segmentation_nii() {
		using PixelType = int;
		using Image = itk::Image<PixelType,3>;
		using Reader = itk::ImageFileReader<Image>;

		Reader::Pointer reader_p = Reader::New();
		reader_p->SetFileName(segmentationPath_nii);
		reader_p->Update();
		Image::Pointer image = reader_p->GetOutput();
		auto spacing = image->GetSpacing();


		//store the metadata dict, so when we have to do ouptut,
		//we may reuse its coordinate system information.
		medDataContainer.medDataDict = image->GetMetaDataDictionary();
		//medDataContainer.medDataDict.Print(std::cout);

		std::cout << "qoffset values read from seg nii file: "<< std::endl;
		if(medDataContainer.medDataDict.HasKey("qoffset_x")) std::cout <<"qoffset_x : "<<medDataContainer.medDataDict.Get("qoffset_x") << std::endl;
		if(medDataContainer.medDataDict.HasKey("qoffset_y")) std::cout <<"qoffset_y : "<<medDataContainer.medDataDict.Get("qoffset_y") << std::endl;
		if(medDataContainer.medDataDict.HasKey("qoffset_z")) std::cout <<"qoffset_z : "<<medDataContainer.medDataDict.Get("qoffset_z") << std::endl;

		medDataContainer.segmentationData_ND_p = std::shared_ptr<SegmentationContainer>(
				new SegmentationContainer(boost::extents[medDataContainer.sizes[0]][medDataContainer.sizes[1]][medDataContainer.sizes[2]]));

		for (unsigned int i = 0; i < medDataContainer.sizes[0]; ++i) {
			for (unsigned int j = 0; j < medDataContainer.sizes[1]; ++j) {
				for (unsigned int k = 0; k < medDataContainer.sizes[2]; ++k) {
					auto pixel = image->GetPixel( { i, j, k });		//get the seg value at the given image position
					((*(medDataContainer.segmentationData_ND_p))[i][j][k]) = pixel;
				}
			}
		}
	}

It seems to store the metadata dictionary, but when I then write my simulation results to file, the offsets and
matrices are off:

void writeNiftiFile(const int &step, const double &t) {
//0. create Buffer,
	unsigned int numberOfPixels = itk_size[0] * itk_size[1];
	if constexpr(DIM ==3 ) numberOfPixels = itk_size[0] * itk_size[1] * itk_size[2];
	//if (master) std::cout << "cellMapper-size: " << cellMapper_p->size() << " and dataset size: " << numberOfPixels << std::endl;
	//if (master) std::cout << "vertex	Mapper-size: " << vertexMapper_p->size() << " and dataset size: " << numberOfPixels << std::endl;
	std::vector<unsigned char> buffer(0);
	buffer.resize(numberOfPixels);

//1. read buffer from simulation data
	//if (master) std::cout << "reading local data into buffer:" << numberOfPixels << std::endl;
	read3DDataIntoBuffer(buffer);
//	int nonzero_entries = 0;
//	for (const auto &pixel : buffer) {
//		if (pixel != 0) nonzero_entries++;
//	}

// each rank sends its buffer to the master rank, which adds
// it to the global one, note that all non-owned fields have
// to be zero!
//auto comm = gv_p->comm();		//class Dune::CollectiveCommunication<ompi_communicator_t*>
	auto comm = mpiHelper_ref.getCollectiveCommunication();    //class Dune::CollectiveCommunication<ompi_communicator_t*>’
	auto * buffer_p = &buffer[0];
	int mpi_status = comm.allreduce<std::plus<double>>(buffer_p, numberOfPixels);
	if (!mpi_status == 0) DUNE_THROW(Dune::Exception, "mpi returned nonzero status in medData output method");

	//2. read buffer into ITK image
	if (master) {
		//std::cout << "nonzero entries after mpi " << ((int) std::count_if(buffer.begin(), buffer.end(), [](auto c) {return c > 0;})) << std::endl;
		//types:
		using PixelType = unsigned char;
		using OutputImageType = itk::Image< PixelType, DIM >;
		using ImportImageType = itk::Image< PixelType, DIM >;
		using ImportFilterType = itk::ImportImageFilter< PixelType, DIM>;
		typename ImportFilterType::Pointer importFilter_p = ImportFilterType::New();

		typename ImportFilterType::IndexType start;
		start.Fill(0);

		typename ImportFilterType::RegionType region;
		region.SetIndex(start);
		region.SetSize(itk_size);
		importFilter_p->SetRegion(region);

		//set origin to middle:
		 itk::SpacePrecisionType itkspacing[DIM] = {0.0};
		for (int i = 0; i < DIM; ++i)
			itkspacing[i] = spacing[i] * 1000;
		importFilter_p->SetSpacing(itkspacing);

		double importOrigin[DIM] = {0.0};
		for (int i = 0; i < DIM; ++i)
			importOrigin[i] = ((double) itk_size[i]) * itkspacing[i] * 0.5;
		importFilter_p->SetOrigin(importOrigin);

//		std::cout << "importFilter spacing after set: " << importFilter_p->GetSpacing() << std::endl;

		const bool importImageFilterWillOwnTheBuffer = false;
		importFilter_p->SetImportPointer(buffer.data(), buffer.size(), importImageFilterWillOwnTheBuffer);
		try {
			importFilter_p->Update();
		} catch (itk::ExceptionObject & e) {
			std::cerr << "exception in file writer " << std::endl;
			std::cerr << e.GetDescription() << std::endl;
			std::cerr << e.GetLocation() << std::endl;
			return;
		}


		auto dict = medData_p->getMedDataContainer().medDataDict;
		std::vector<std::string> keys = dict.GetKeys();
		if(dict.HasKey("qoffset_x")) std::cout <<"qoffset_x : "<<dict.Get("qoffset_x") << std::endl;
		if(dict.HasKey("qoffset_y")) std::cout <<"qoffset_y : "<<dict.Get("qoffset_y") << std::endl;
		if(dict.HasKey("qoffset_z")) std::cout <<"qoffset_z : "<<dict.Get("qoffset_z") << std::endl;

		importFilter_p->SetMetaDataDictionary(medData_p->getMedDataContainer().medDataDict);
		typedef itk::ImageFileWriter<OutputImageType> WriterType;

		//prepare folder
		std::ostringstream newFolderName;
		newFolderName << "step" << std::setw(5) << std::setfill('0') << step;
		std::string subFolder = createFolder(pTree.get < std::string > ("output.medical_format_output_path"), newFolderName.str());  //
		auto filename = subFolder + "/" + "output3D.nii.gz";

		//write file
		typename WriterType::Pointer writer = WriterType::New();
		writer->SetFileName(filename);
		writer->SetInput(importFilter_p->GetOutput());
		writer->Update();
	}
}

Is there some fancy hidden transformation that is beeing applied? Is there documentation about this behaviour?

Any Help is greatly appreciated!

Michael

If you want to preserve all metadata in one fell swoop, use CopyInformation, e.g. outputImage->CopyInformation(originalImage);

In your code, you seem to be forgetting to copy direction, e.g. importFilter_p->SetDirection(image->GetDirection());

1 Like