Well, as is often obnoxiously said in conference talks – this is not really a question, more of a comment…
One of the things that raise worries when working with ITK from Python, is the issue of types. In most cases, Python container types (lists, tuples etc) get translated to the proper C++ types (vectors, matrices etc) automatically, but there are edge cases, and it feels safer to get an object from the library and modify it, rather than try to create a new object in Python and hope it fits.
I was trying to turn a set of images (PNGs for example) to DICOM files; read the series into a volume, then write the volume as a series. This mostly involves setting up the proper metadata for DICOM. But there is also something about the images themselves: They do not really have geometric information, just number of pixels. In order for the read volume to make sense – even assuming the images are all aligned and ordered – we need to set the spacing.
In the spirit of the first paragraph above, I did this like so:
reader = itk.ImageSeriesReader[ImageType].New()
# ... set up reader
reader.Update() # Actual reading
image_3d = reader.GetOutput()
spacing = image_3d.GetSpacing()
spacing[0] = spacing[1] = pixel_size
spacing[2] = z_space_between_images
image_3d.SetSpacing(spacing)
And this looked like it was working, but it really wasn’t – the x and y spacing worked, but the z spacing behaved as if it was never set. When DICOM images were written, their z-origin progressed in steps of 1.0, and not the value I had selected.
Only after I dug in, and played more with it, I found the reason. Without getting into too much details, correct z stepping relies on some internal calculation that is performed when the spacing is changed; but ITK does not perform it when SetSpacing() is called with the value it already has, because that would be redundant.
Except that, the above “try to be safe” method makes it so that the spacing vector is changed in place – which means, that when I later call SetSpacing(), it is indeed with the same vector as it already has. Optimization kicks in, the calculation is skipped, z spacing is broken.
The fix:
reader = itk.ImageSeriesReader[ImageType].New()
# ... set up reader
reader.Update() # Actual reading
image_3d = reader.GetOutput()
spacing = [pixel_size, pixel_size, z_space]
image_3d.SetSpacing(spacing)
I’m sure this case is not so special, but common to many places in ITK. Thought this might save some time for some other poor souls.