ImageData.zoom_image() size orders

AnderBiguri commented 4 years ago

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 commented 4 years ago

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

AnderBiguri commented 4 years ago

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

KrisThielemans commented 3 years ago

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

AnderBiguri commented 3 years ago

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

AnderBiguri commented 3 years ago

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 commented 3 years ago

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 commented 3 years ago

@KrisThielemans correct.

KrisThielemans commented 3 years ago

I'm very confused after looking at https://github.com/SyneRBI/SIRF/blob/b36c1b2378644d6cbc2b4cdab53210e3e16db25c/src/common/GeometricalInfo.cpp#L43-L49 as this seems alright (noting that offset etc are in LPS).

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

DANAJK commented 3 years ago

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 commented 3 years ago

Strange - can you quickly do:

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

ashgillman commented 3 years ago

I'm looking at https://github.com/SyneRBI/SIRF/blob/2cd7275dcc0c09a01e30314cf95a7065405c3aa7/src/common/SIRF.py#L683-L699

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: https://github.com/SyneRBI/SIRF/blob/720ed3d55d25f46846060746d3ba41a353ed1ed6/src/xSTIR/cSTIR/stir_data_containers.cpp#L589-L624

ashgillman commented 3 years ago

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 https://github.com/SyneRBI/SIRF/commit/bb7690e0c333d673c4f3c712696e45ea827b542b 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 commented 3 years ago

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 commented 3 years ago

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?

SyneRBI / SIRF

ImageData.zoom_image() size orders #812