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