Orientation in Diffusion Tensor Imaging file (DTI) .nii


#1

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


(Dženan Zukić) #2

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());