ImageData.zoom_image() size orders

[AnderBiguri]

I am trying to change the geometric info of an sirf.STIR.ImageData object.

While doing some tests, I tried:

import sirf.STIR as stir
stir.ImageData("myimage.hv")

off=img.get_geometrical_info().get_offset()

img.get_geometrical_info().print_info()
img=img.zoom_image(offsets_in_mm=(0.0,0.0,-off[2]))
img.get_geometrical_info().print_info()

Which prints:

Offset: (-350, -350, 21.52)
Spacing: (1.36719, 1.36719, 2.80005)
Size: (512, 512, 411)
Dir mat: 
1, 0, 0
0, 1, 0
0, 0, -1

Offset: (-371.52, -350, 21.52)
Spacing: (1.36719, 1.36719, 2.80005)
Size: (512, 512, 411)
Dir mat: 
1, 0, 0
0, 1, 0
0, 0, -1

What threw me off is that offset is printed as (-350, -350, 21.52), the variable off has values off=(-350, -350, 21.52), yet, when I use the last value of off into offsets_in_mm, its the first value of the offset that gets changed.

This seems just a mismatch of python vs other (C++?) dimension storage issue (ImageData is indexed (z,x,y)). But it creates a strong undocumented inconsistency. I would expect that if I get the offset from an object, then set it to its negative, to return zero, without the need of rearranging the offset variable order.

Note this happens with all variables, not only offsets_in_mm.

I guess this is either the same, or a effect of #791

[KrisThielemans]

duplicate or similar to #791. we do need to sort this out of course.

[AnderBiguri]

@KrisThielemans yes indeed, I already linked it in the issue. Did not realize it was a similar thing until after I posted it.

[KrisThielemans]

@AnderBiguri could you have a look at this to see if this is required for v3.0? I can't remember...

[AnderBiguri]

I think this is closed with #791, let me quick test.

[AnderBiguri]

So this is the current stage of this:

off=image.get_geometrical_info().get_offset()
print(off)

image.get_geometrical_info().print_info()
image=image.zoom_image(offsets_in_mm=(0.0,0.0,-off[2]))
image.get_geometrical_info().print_info()

returns

(-0.0, -407.4256, -407.4256)

Offset: (-407.426, -407.426, -0)
Spacing: (2.7344, 2.7344, 2.76148)
Size: (298, 298, 660)
Dir mat: 
1, 0, 0
0, 1, 0
0, 0, -1


Offset: (0, -407.426, -0)
Spacing: (2.7344, 2.7344, 2.76148)
Size: (298, 298, 660)
Dir mat: 
1, 0, 0
0, 1, 0
0, 0, -1

I still think this is confusing because of the printing order, but at least the indices correspond to the correct ones, thanks to #791

[KrisThielemans]

thanks. So is the remaining problem then that print_info() prints the offsets (and I guess spacing and size) in the wrong order (for C++ and Python)?

[AnderBiguri]

@KrisThielemans correct.

[KrisThielemans]

I'm very confused after looking at

SIRF/src/common/GeometricalInfo.cpp

Lines 43 to 49 in b36c1b2

	info << _offset[0] << ", " << _offset[1] << ", " << _offset[2] << ")\n";
	info << "Spacing: (";
	info << _spacing[0] << ", " << _spacing[1] << ", " << _spacing[2] << ")\n";
	info << "Size: (";
	info << _size[0] << ", " << _size[1] << ", " << _size[2] << ")\n";

as this seems alright (noting that offset etc are in LPS).

@evgueni-ovtchinnikov @ashgillman do you understand what's going on?

[DANAJK]

offset is the location in LPS space of the voxel that has image coordinate (0,0,0). offset is in LPS space so the ordering should never change.

The word 'zoom' implies 'a closer look' so I am not sure what zoom_image is doing, or why it takes offset as an input.

PixelSpacing in DICOM is a 2-element vector with meaning [height width] in mm for a pixel. I think the suggestion above is that the order in spacing should not be tied to height or width, but apply to the corresponding array dimension (this is like NIfTI). For example a spacing of [ 4 5 6 ] would mean that for an array indexed as img(i1, i2, i3), then if you increase the value of i2 by 1, you move 5mm in the direction of increasing i2. Incidentally, this direction is the vector which is the 2nd column (counting from 1) of the dir matrix. The elements of this column vector are always in the order LPS.

So, if the array dimensions are swapped around:

• offset should not change - it is in LPS and points to (0,0,0)
• the elements of spacing and the columns of the dir matrix should correspond in order to the indices into the array.
• the rows of the dir matrix are LPS and should not be swapped.
• it needs to be clear if the reported parameters are referring to an internal SIRF representation, or the as_array returned in Python or MATLAB.
• for completeness, the mapping between data (conceptually a 1D line of data in a file or a stream) and the array needs to be specified. The NIfTI way would be that data is placed into the array with the first index varying most rapidly.

[ashgillman]

Strange - can you quickly do:

with open("myimage.hv", 'r') as fp:
     print(fp.read())

[ashgillman]

I'm looking at

SIRF/src/common/SIRF.py

Lines 683 to 699 in 2cd7275

	def get_offset(self):
	"""Offset is the LPS coordinate of the centre of the first voxel."""
	arr = numpy.ndarray((3,), dtype = numpy.float32)
	try_calling(pysirf.cSIRF_GeomInfo_get_offset(self.handle, arr.ctypes.data))
	return tuple(arr[::-1])
	def get_spacing(self):
	"""Spacing is the physical distance between voxels in each dimension."""
	arr = numpy.ndarray((3,), dtype = numpy.float32)
	try_calling (pysirf.cSIRF_GeomInfo_get_spacing(self.handle, arr.ctypes.data))
	return tuple(arr[::-1])
	def get_size(self):
	"""Size is the number of voxels in each dimension."""
	arr = numpy.ndarray((3,), dtype = numpy.int32)
	try_calling (pysirf.cSIRF_GeomInfo_get_size(self.handle, arr.ctypes.data))
	return tuple(arr[::-1])

I'm not sure originally why these weren't returning LPS, but I think this is just obfuscating it. We want GeometricalInfo to be always defined and stored in LPS and it should be defined that way when creating the ImageData from a STIR object.

Actually, I thought we did do that, and looking at the code, we certainly tried:

SIRF/src/xSTIR/cSTIR/stir_data_containers.cpp

Lines 589 to 624 in 720ed3d

	// SIRF offest is STIR's LPS location of first voxel
	VoxelisedGeometricalInfo3D::Offset offset;
	const stir::CartesianCoordinate3D<int> first_vox
	= vox_image->get_min_indices();
	const Coord3DF first_vox_coord
	= vox_image->get_LPS_coordinates_for_indices(first_vox);
	offset[0] = first_vox_coord.x();
	offset[1] = first_vox_coord.y();
	offset[2] = first_vox_coord.z();
	// SIRF and STIR share size definition
	VoxelisedGeometricalInfo3D::Size size;
	size[0] = vox_image->get_x_size();
	size[1] = vox_image->get_y_size();
	size[2] = vox_image->get_z_size();
	// SIRF's spacing is STIR's voxel size, but with different order
	VoxelisedGeometricalInfo3D::Spacing spacing;
	spacing[0] = vox_image->get_voxel_size()[3];
	spacing[1] = vox_image->get_voxel_size()[2];
	spacing[2] = vox_image->get_voxel_size()[1];
	// Find axes direction as the normalised vector between voxels
	// in each direction
	VoxelisedGeometricalInfo3D::DirectionMatrix direction;
	for (int axis = 0; axis < 3; axis++) {
	Coord3DI next_vox_along_axis(first_vox);
	next_vox_along_axis[3 - axis] += 1;
	const Coord3DF next_vox_along_axis_coord
	= vox_image->get_LPS_coordinates_for_indices(next_vox_along_axis);
	Coord3DF axis_direction
	= next_vox_along_axis_coord - first_vox_coord;
	axis_direction /= stir::norm(axis_direction);
	for (int dim = 0; dim < 3; dim++)
	direction[dim][axis] = axis_direction[3 - dim];
	}

[ashgillman]

I'm a bit confused by #791 actually, I am currently using SIRF v2.2.0, which is older than that issue, and my geometrical info orders are all correct /except/ for the zoom_images as outlined in this issue. There is no MWE in that issue - I think that the issue here is caused by bb7690e which be reverted.

template_PET_raw_path = os.path.join(
    examples_data_path('PET'), 'mMR/mMR_template_span11.hs')
template_PET_raw = pet.AcquisitionData(template_PET_raw_path)
template_PET_im = pet.ImageData(template_PET_raw)
print('Template:')
print('  shape:', template_PET_im.get_geometrical_info().get_size())
print('  spacing:', template_PET_im.get_geometrical_info().get_spacing())
print('  origin:', template_PET_im.get_geometrical_info().get_offset())

Template:
  shape: (285, 285, 127)
  spacing: (2.08626, 2.08626, 2.03125)
  origin: (-296.24893, -296.24893, -0.0)

[AnderBiguri]

Strange - can you quickly do:

with open("myimage.hv", 'r') as fp:
     print(fp.read())

!INTERFILE :=
originating system := Discovery MI
!imaging modality := PT
!version of keys := STIR4.0
name of data file := dynamicXcat_WB_1260_1295.v
!GENERAL DATA :=
!GENERAL IMAGE DATA :=
!type of data := PET
imagedata byte order := LITTLEENDIAN
!PET STUDY (General) :=
!PET data type := Image
!number format := float
!number of bytes per pixel := 4
number of dimensions := 3
matrix axis label [1] := x
!matrix size [1] := 298
scaling factor (mm/pixel) [1] := 2.7344
matrix axis label [2] := y
!matrix size [2] := 298
scaling factor (mm/pixel) [2] := 2.7344
matrix axis label [3] := z
!matrix size [3] := 660
scaling factor (mm/pixel) [3] := 2.76148
first pixel offset (mm) [1] := -407.4256
first pixel offset (mm) [2] := -407.4256
first pixel offset (mm) [3] := 0
number of time frames := 1
image duration (sec)[1] := 35
image relative start time (sec)[1] := 1260
number of energy windows := 1
energy window lower level[1] := 650
energy window upper level[1] := 425
!END OF INTERFILE :=

[AnderBiguri]

For #791 testing, can you try this (on your prior version):

import sirf.STIR as pet
acd=pet.AcquisitionData("Discovery MI")
image=acd.create_uniform_image()
off=image.get_geometrical_info().get_offset()
print(off)
image.get_geometrical_info().print_info()
image=image.zoom_image(offsets_in_mm=(0.0,0.0,-off[2]))
image.get_geometrical_info().print_info()

Now, it outputs:

(-0.0, -355.16602, -355.16602)

Offset: (-355.166, -355.166, -0)
Spacing: (2.206, 2.206, 2.76148)
Size: (323, 323, 53)
Dir mat: 
1, 0, 0
0, 1, 0
0, 0, -1


Offset: (0, -355.166, -0)
Spacing: (2.206, 2.206, 2.76148)
Size: (323, 323, 53)
Dir mat: 
1, 0, 0
0, 1, 0
0, 0, -1

I believe that the problem was that off would the be ( -355.16602, -355.16602, -0.0), so if you wanted to change that first offset (i.e. replicate what the code above does), you'd need to do:

image=image.zoom_image(offsets_in_mm=(0.0,0.0,-off[0]))

Which is confusing, as you need to pass index [0] to change the first dim, as the third dim to zoom_image

The problem may be that we solved #791 wrong. Instead of changing what get_offset() returned, maybe we needed to change the order of the inputs in zoom_image?