[Insight-users] ExtractImageFilter & ChangeInformationImageFilter to keep slice origin

[Antoine Leimgruber]

Hello,

I have a 3D image and I am writing it slice by slice in the z direction using ExtractImageFilter. I use an ImageFileWriter and call my files 0001.dcm, 0002.dcm,....

What is the best way of not loosing the origin of the slice when writing it to 000i.dcm ?

This is what I tried

I first used ExtractImageFilter as in the ITK examples, but I noticed that every extracted slice had the origin [0 0] and also that there was no image position patient tag in the dicom file.

To have a 3D origin, I defined my "slice" as a Dimension=3 image and then I set the size[2] of the region to extract to 1. With that, every extracted slice had the origin [0 0 0].

I used then a ChangeInformationFilter to change the origin of the extracted slice. With the spacing and the slice index I thought that it was an easy solution.

However, at the end I do not see a "ImagePositionPatient" (0020,0032) tag in the header of my slices, so there is no origin at all.

Here is the code,

Thank you for your help

Antoine


/*
-----------------------------------------------------------------
	write dicom series from my volume
-----------------------------------------------------------------
*/	
	// declare a new Image type for the slice
	const unsigned int	SliceDimension = 3; // changed 2 to 3 to have a 3D origin
	typedef itk::Image< unsigned short, SliceDimension > SliceType;
	SliceType::Pointer slice = SliceType::New();

	// attempt to manually modify the origin
	typedef itk::ChangeInformationImageFilter<ImageType> ChangeFilterType4Slice;
	ChangeFilterType4Slice::Pointer sliceheaderfilter = ChangeFilterType4Slice::New();
	double neworigin [SliceDimension];
	neworigin[0]=0;
	neworigin[1]=0;
	neworigin[2]=0; // will be updated in loop for each slice

	// declare a new writer type for writing slices
	typedef itk::ImageFileWriter< SliceType > SliceWriterType;
	SliceWriterType::Pointer slicewriter = SliceWriterType::New();
	
	// declare an extract filter to extract slices
	typedef itk::ExtractImageFilter< ImageType, SliceType > ExtractSliceFilterType;
	ExtractSliceFilterType::Pointer  extractslicefilter = ExtractSliceFilterType::New();
	
	// connect the input of the filter to the output of the previous filter
	if (fR) {extractslicefilter->SetInput( samplfilter->GetOutput() );}
	else{extractslicefilter->SetInput( headerfilter->GetOutput() );}
	
	// define a region to extract from data coming out of previous filters
	ImageType::RegionType extractregion;
	if (fR) {extractregion = samplfilter->GetOutput()->GetLargestPossibleRegion();}
	else{extractregion = headerfilter->GetOutput()->GetLargestPossibleRegion();}
	
	// define that we want to extract a slice in the k direction
    ImageType::SizeType extractsize = extractregion.GetSize();
	extractsize[2] = 1; // I changed that when changing SliceDimension from 2 to 3
	
	// define the size of the slice
	ImageType::RegionType desiredRegion;
	desiredRegion.SetSize(  extractsize  );
	
	// set a variable indexing which slice will be extracted
	ImageType::IndexType extractstart;
	
	// initialize the indexing variable
	extractstart[0]=0; 
	extractstart[1]=0; 
	extractstart[2]=0;
	
	// create a name for the slices
	std::string SliceName;
	int SliceNumber;
	
	// when the loop will end
	const int smax = (extractregion.GetSize())[2];
	
	// for each slice
	for (int s = 0; s < smax; s++) 
	{
		
		// extract slice
		extractstart[2]=s; // redundant for k=0
		desiredRegion.SetIndex( extractstart );
		std::cout << "Extract     " <<extractstart << " into file ";
		extractslicefilter->SetExtractionRegion( desiredRegion );

		// build a file name (quick fix)
		char myBuf[20];
		SliceNumber=s+1;
		sprintf(myBuf, "%i", SliceNumber); 
		SliceName = "";
		if (SliceNumber < 10) {SliceName += "000";}
		else if (SliceNumber <100) {SliceName += "00";}
		else if (SliceNumber <1000) {SliceName += "0";}
		else if (SliceNumber <10000) {SliceName += "";}
		else {
				std::cout << "Too many slices !\n";
				return EXIT_FAILURE;
		}
		SliceName += myBuf;
		SliceName += ".dcm";
		std::cout << SliceName;
		
		
		slicewriter->SetFileName( SliceName );

		// change origin changeinformationfilter 
		sliceheaderfilter->SetInput( extractslicefilter->GetOutput() );
		neworigin[2]=s*spacing[2];
		sliceheaderfilter->SetOutputOrigin(neworigin);
		sliceheaderfilter->ChangeOriginOn();
		sliceheaderfilter->Update();
		slicewriter->SetInput( sliceheaderfilter->GetOutput() );
		

		try 
		{ 
		slicewriter->Update(); 
		} 
		catch( itk::ExceptionObject & err ) 
		{ 
		std::cerr << "ExceptionObject caught !" << std::endl; 
		std::cerr << err << std::endl; 
		return EXIT_FAILURE;
		} 
		std::cout << " with origin " << sliceheaderfilter->GetOutput()->GetOrigin() << std::endl;
		
	}