Insert a 2D slice in a 3D image

Hi all,
I created a 3D image with desirable spacing and origin, now I want to insert a 2D slice in one of the slices of this 3D image. Does anyone know how I can do it? Below, I also sent a code written for creating 3D image.

// Software Guide : BeginCodeSnippet
#include "itkImage.h"
// Software Guide : EndCodeSnippet
#include "itkImageFileWriter.h"
int
main(int, char *[])
{
  // Software Guide : BeginLatex
  //
  // Then we must decide with what type to represent the pixels
  // and what the dimension of the image will be. With these two
  // parameters we can instantiate the \code{Image} class. Here we create
  // a 3D image with \code{unsigned short} pixel data.
  //
  // Software Guide : EndLatex
  //
  // Software Guide : BeginCodeSnippet
  using ImageType = itk::Image<unsigned short, 3>;
  // Software Guide : EndCodeSnippet
  using WriterType = itk::ImageFileWriter<ImageType>;
 
  // Software Guide : BeginLatex
  //
  // The image can then be created by invoking the \code{New()} operator
  // from the corresponding image type and assigning the result
  // to a \doxygen{SmartPointer}.
  //
  // \index{itk::Image!Pointer}
  // \index{itk::Image!New()}
  //
  // Software Guide : EndLatex
  //
  // Software Guide : BeginCodeSnippet
  ImageType::Pointer image = ImageType::New();
  // Software Guide : EndCodeSnippet
 
 
  // Software Guide : BeginLatex
  //
  // In ITK, images exist in combination with one or more
  // \emph{regions}. A region is a subset of the image and indicates a
  // portion of the image that may be processed by other classes in
  // the system. One of the most common regions is the
  // \emph{LargestPossibleRegion}, which defines the image in its
  // entirety. Other important regions found in ITK are the
  // \emph{BufferedRegion}, which is the portion of the image actually
  // maintained in memory, and the \emph{RequestedRegion}, which is
  // the region requested by a filter or other class when operating on
  // the image.
  //
  // In ITK, manually creating an image requires that the image is
  // instantiated as previously shown, and that regions describing the image are
  // then associated with it.
  //
  // A region is defined by two classes: the \doxygen{Index} and
  // \doxygen{Size} classes. The origin of the region within the
  // image is defined by the \code{Index}. The extent, or size, of the region
  // is defined by the \code{Size}. When an image is created manually, the user
  // is responsible for defining the image size and the index at which the
  // image grid starts. These two parameters make it possible to process
  // selected regions.
  //
  // The \code{Index} is represented by a n-dimensional array where each
  // component is an integer
  // indicating---in topological image coordinates---the initial pixel
  // of the image.
  //
  // \index{itk::Image!Size}
  // \index{itk::Image!SizeType}
  //
  // Software Guide : EndLatex
  //
  // Software Guide : BeginCodeSnippet
  ImageType::IndexType start;
 
  start[0] = 0; // first index on X
  start[1] = 0; // first index on Y
  start[2] = 0; // first index on Z
  // Software Guide : EndCodeSnippet
 
  // Software Guide : BeginLatex
  //
  // The region size is represented by an array of the same dimension as the
  // image (using the \doxygen{Size} class). The components of the array are
  // unsigned integers indicating the extent in pixels of the image along
  // every dimension.
  //
  // \index{itk::Image!Index}
  // \index{itk::Image!IndexType}
  //
  // Software Guide : EndLatex
  //
  // Software Guide : BeginCodeSnippet
  ImageType::SizeType size;
 
  size[0] = 148; // size along X
  size[1] = 224; // size along Y
  size[2] = 42; // size along Z
  // Software Guide : EndCodeSnippet
 
  // Software Guide : BeginLatex
  //
  // Having defined the starting index and the image size, these two
  // parameters are used to create an \doxygen{ImageRegion} object
  // which basically encapsulates both concepts.
  // The region is initialized with the starting index and size of the
  // image.
  //
  // \index{itk::Image!itk::ImageRegion}
  // \index{itk::Image!RegionType}
  //
  // Software Guide : EndLatex
  double spacing[3];
  spacing[0] = 1.00390625;
  spacing[1] = 1.00390625;
  spacing[2] = 3.00000000;

  double origin[3];
  origin[0] = -155.60546875;
  origin[1] = -29.42187500;
  origin[2] = -648.00000000;
 
  // Software Guide : BeginCodeSnippet
  ImageType::RegionType region;
 
  region.SetSize(size);
  region.SetIndex(start);
  
  // Software Guide : EndCodeSnippet
 
  // Software Guide : BeginLatex
  //
  // Finally, the region is passed to the \code{Image} object in order to define its
  // extent and origin. The \code{SetRegions} method sets the
  // \emph{LargestPossibleRegion}, \emph{BufferedRegion}, and \emph{RequestedRegion}
  // simultaneously. Note that none of the operations performed to this point
  // have allocated memory for the image pixel data. It is necessary to
  // invoke the \code{Allocate()} method to do this. Allocate does not
  // require any arguments since all the information needed for memory
  // allocation has already been provided by the region.
  //
  // \index{itk::Image!Allocate()}
  // \index{itk::Image!SetRegions()}
  //
  // Software Guide : EndLatex
 
  // Software Guide : BeginCodeSnippet
  image->SetRegions(region);
  
  image->Allocate();
  // Software Guide : EndCodeSnippet
  image->SetSpacing(spacing);
  image->SetOrigin(origin);

  WriterType::Pointer writer = WriterType::New();
  writer->SetFileName("output.nrrd");
  writer->SetInput(image);
  try
  {
	  writer->Update();
  }
  catch (const itk::ExceptionObject & err)
  {
	  std::cerr << "ExceptionObject caught !" << std::endl;
	  std::cerr << "err" << std::endl;
	  return EXIT_FAILURE;
  }
 
  return EXIT_SUCCESS;
}

The normal way would be using iterators. You would need a loop similar to this one, but set up to iterate over your slice, and desired slice in the 3D image. Additional examples: 1, 2.

Alternatively, you could find out location of memory belonging to desired slice, and do a memcpy or std::copy. Key method: GetBufferPointer().

2 Likes

I answered a similar question on StackOverflow a while back:

I did it in SimpleITK/Python by converting my 2d slice into a 3d volume (with zdim=1) and then using the Paste function. You can do the same thing in ITK with the PasteImageFilter:

https://itk.org/Doxygen/html/classitk_1_1PasteImageFilter.html

1 Like

Hi,
Thanks for your reply. What I did, I first create a 3D image with desirable origin and spacing, then based on your suggestion, I converted my 2D image to 3D image with dimz=1, and I wrote the following code to paste that image to desirable slice( first image) but now, I am just seeing a white image in the first slice instead of desirable slice. Can you please help me with the code and see if there is a problem with that?

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkPasteImageFilter.h"

int
main(int argc, char * argv[])
{
	if (argc != 3)
	{
		std::cerr << "Usage: " << std::endl;
		std::cerr << argv[0];
		std::cerr << " <SourceFileName> <DestinationFileName> <OutputFileName>";
		std::cerr << std::endl;
		return EXIT_FAILURE;
	}

	const char * sourceFileName = argv[1];
	const char * destinationFileName = argv[2];
	const char * outputFileName = argv[3];

	/*int startX = std::stoi(argv[4]);
	int startY = std::stoi(argv[5]);*/

	constexpr unsigned int Dimension = 3;

	using PixelType = unsigned short;
	using ImageType = itk::Image<PixelType, Dimension>;

	ImageType::IndexType index;

	index[2] = 1;

	using ReaderType = itk::ImageFileReader<ImageType>;
	ReaderType::Pointer sourceReader = ReaderType::New();
	sourceReader->SetFileName(sourceFileName);
	sourceReader->Update();

	ReaderType::Pointer destinationReader = ReaderType::New();
	destinationReader->SetFileName(destinationFileName);

	using FilterType = itk::PasteImageFilter<ImageType, ImageType>;
	FilterType::Pointer filter = FilterType::New();
	filter->SetSourceImage(sourceReader->GetOutput());
	filter->SetSourceRegion(sourceReader->GetOutput()->GetLargestPossibleRegion());
	filter->SetDestinationImage(destinationReader->GetOutput());
	filter->SetDestinationIndex(index);

	using WriterType = itk::ImageFileWriter<ImageType>;
	WriterType::Pointer writer = WriterType::New();
	writer->SetFileName(outputFileName);
	writer->SetInput(filter->GetOutput());
	try
	{
		writer->Update();
	}
	catch (itk::ExceptionObject & error)
	{
		std::cerr << "Error: " << error << std::endl;
		return EXIT_FAILURE;
	}

	return EXIT_SUCCESS;
}

You didn’t initialize the first two elements of index. So when I ran your code they were getting random values.

I added the following lines but there was not effect.

index[0] = 0;
index[1] = 0;

Do you have another suggestion?

I don’t know. It worked for me. I took your code, put in those two lines, and tried it with a test slice and volume.

BTW, the first if should be if (argc != 4)

Is it possible for you to test the code on the attached images. I want to see the source image in one of the slices of destination image. I really appreciate it. images.rar (2.7 KB)

It works. But visually seeing the difference is nearly impossible because your source slice has values 0 and 1, so its intensities are completely swamped by 50000 or so of the destination pixels.

#include "itkImageFileReader.h"
#include "itkImageFileWriter.h"
#include "itkPasteImageFilter.h"

int
main(int argc, char * argv[])
{
  if (argc < 4)
  {
    std::cerr << "Usage: " << std::endl;
    std::cerr << argv[0];
    std::cerr << " <SourceFileName> <DestinationFileName> <OutputFileName>";
    std::cerr << std::endl;
    return EXIT_FAILURE;
  }

  const char * sourceFileName = argv[1];
  const char * destinationFileName = argv[2];
  const char * outputFileName = argv[3];

  /*int startX = std::stoi(argv[4]);
  int startY = std::stoi(argv[5]);*/

  constexpr unsigned int Dimension = 3;

  using PixelType = unsigned short;
  using ImageType = itk::Image<PixelType, Dimension>;

  ImageType::IndexType index;
  index[0] = 0;
  index[1] = 0;
  index[2] = 1;

  using ReaderType = itk::ImageFileReader<ImageType>;
  ReaderType::Pointer sourceReader = ReaderType::New();
  sourceReader->SetFileName(sourceFileName);
  sourceReader->Update();

  ReaderType::Pointer destinationReader = ReaderType::New();
  destinationReader->SetFileName(destinationFileName);

  using FilterType = itk::PasteImageFilter<ImageType, ImageType>;
  FilterType::Pointer filter = FilterType::New();
  filter->SetSourceImage(sourceReader->GetOutput());
  filter->SetSourceRegion(sourceReader->GetOutput()->GetLargestPossibleRegion());
  filter->SetDestinationImage(destinationReader->GetOutput());
  filter->SetDestinationIndex(index);

  using WriterType = itk::ImageFileWriter<ImageType>;
  WriterType::Pointer writer = WriterType::New();
  writer->SetFileName(outputFileName);
  writer->SetInput(filter->GetOutput());
  try
  {
    writer->Update();
  }
  catch (itk::ExceptionObject & error)
  {
    std::cerr << "Error: " << error << std::endl;
    return EXIT_FAILURE;
  }

  return EXIT_SUCCESS;
}

Yeah, the source image is a binary mask (1’s and 0’s), and the destination image pixels all have the same value, 52685.