get data container from acquisition model

mehrhardt commented 6 years ago

I want to get a data container from the acquisition model (see below) which causes an error in STIR. Any idea what is going wrong here?

import os
import shutil
import pSTIR as pet
os.chdir(pet.petmr_data_path('pet'))
shutil.rmtree('working_folder/thorax_single_slice',True)
shutil.copytree('thorax_single_slice','working_folder/thorax_single_slice')
os.chdir('working_folder/thorax_single_slice')

image = pet.ImageData('emission.hv');
image_array=image.as_array()*.05
image.fill(image_array);
am = pet.AcquisitionModelUsingRayTracingMatrix()
am.set_num_tangential_LORs(5)
templ = pet.AcquisitionData('template_sinogram.hs')
am.set_up(templ,image)

y = am.acq_templ.copy()

File: /home/sirfuser/devel/install/python/sirf/STIR.py
Line: 791
check_status found the following message sent from the engine:
Traceback (most recent call last):

  File "<ipython-input-5-c612a74e83a1>", line 1, in <module>
    y = am.acq_templ.copy()

  File "/home/sirfuser/devel/install/python/sirf/STIR.py", line 333, in copy
    return self.clone()

  File "/home/sirfuser/devel/install/python/sirf/STIR.py", line 811, in clone
    ad.fill(self)

  File "/home/sirfuser/devel/install/python/sirf/STIR.py", line 791, in fill
    (self.handle, value.handle))

  File "/home/sirfuser/devel/install/python/sirf/Utilities.py", line 321, in try_calling
    check_status(returned_handle, inspect.stack()[1])

  File "/home/sirfuser/devel/install/python/sirf/Utilities.py", line 317, in check_status
    raise error(errorMsg)

error: ??? "'ProjDataFromStream: error reading data' exception caught at line 578 of /home/sirfuser/devel/build/sources/SIRF/src/xSTIR/cSTIR/cstir.cpp\nthe reconstruction engine output may provide more information"

paskino commented 6 years ago

Apparently you are trying to create an AcquisitionData from a template (and probably fill it with zeros?). The right way to do should be

y = pet.AcquisitionData(am.acq_templ)
y.fill(0)

What goes wrong in y = am.acq_templ.copy() is that you are copying a template which doesn't contain data: error: ??? "'ProjDataFromStream: error reading data'. Although a template inherits from AcquisitionData it is meant to contain only metadata not data.

mehrhardt commented 6 years ago

Yes, that makes perfect sense. I will use it.

However, it is still strange that am.acq_templ has a copy function but one should not use it ...

paskino commented 6 years ago

You are using the AcquisitionModel (Operator in CCPi terms) to create an instance of AcquisitionData (or ImageData) it operates between.

In CIL we defined 2 ("virtual") methods in the Operator to do this, allocate_direct and allocate_adjoint.

However, we discussed with @evgueni-ovtchinnikov about this issue and his opinion is that it should be the job of the acquisition or image template (AcquisitionGeometry or ImageGeometry for the CIL) to allocate an empty instance of ImageData or AcquisitionData, see https://github.com/vais-ral/CCPi-Framework/issues/149. Maybe along these lines,

y = am.acq_templ.allocate_new(0) # the 0 fills the new instance with zeros

In SIRF the template is derived from DataContainer which has the copy method. One should not use the method copy because the template doesn't contain data by definition.

KrisThielemans commented 6 years ago

a few comments.

I thought we already have a create_empty_copy function or whatever the name is (as opposed to allocate_new).
I think that it is reasonable for copy to fail on an empty AcquisitionData, but I suppose it is not obvious to the user that the supplied template AcquisitionData will be empty. I guess we could let STIR throw a more reasonable exception (check that there is no data in there, then write something). would that help?
I am surprised that AcquisitionModel exposes a member acq_templ. I guess this is because Python doesn't really have private variables, but @evgueni-ovtchinnikov, please rename this to __acq_template such that users have to work hard to see it.

mehrhardt commented 6 years ago

I could not find a create_empty_copy function. For the ImageData I found am.img_templ.get_uniform_copy() which is almost doing what I want but unnecessarily filling it with 1s (which the function name does not indicate). In contast, am.acq_templ.create_uniform_image() creates an image of 0s. Finally, for AcquisitionData the same function exists: a = am.acq_templ.get_uniform_copy(), however, now creating a copy with 0s.

Maybe the names should be changed to create_image_ones and create_copy_zeros?

Interestingly, along the lines of what @paskino proposed for AcquisitionData I tried x = pet.ImageData(am.img_templ) but this throws an error:

Traceback (most recent call last):

  File "<ipython-input-27-57190e3669d4>", line 1, in <module>
    x = pet.ImageData(am.img_templ)

  File "/home/sirfuser/devel/install/python/sirf/STIR.py", line 364, in __init__
    ('wrong argument ' + repr(arg) + ' for ImageData constructor')

error: ??? 'wrong argument <sirf.STIR.ImageData object at 0x7f48260dd190> for ImageData constructor'

paskino commented 6 years ago

Does the equivalent for the AcquisitionData

y = pet.AcquisitionData(am.acq_templ)

work?

mehrhardt commented 6 years ago

Yes, this works well.

In the mean time, I tried to use am.acq_templ.create_uniform_image() which failed as it creates an image but not of the dimensions needed. Thinking about this again, why would the AcquisitionData know what image I want? Why does this function even exists in the first place?

KrisThielemans commented 6 years ago

I believe get_uniform_copy has an optional parameter with the value, @evgueni-ovtchinnikov, is this correct ? (should be in the doc.)

mehrhardt commented 6 years ago

I believe get_uniform_copy has an optional parameter with the value, @evgueni-ovtchinnikov, is this correct ? (should be in the doc.)

Is there an easy way to see the doc in spyder?

KrisThielemans commented 6 years ago

RTM :-)

In spyder, ? somemethod should work.

KrisThielemans commented 6 years ago

why would the AcquisitionData know what image I want? Why does this function even exists in the first place?

PET images are determined by what was scanned. Currently, this gives for instance z-spacing and number of slices (but also transaxial fov). Hopefully soon this will also take bed position into account.

In MR, this is even more the case, as with FFT, you have no choice about voxel size.

paskino commented 5 years ago

@evgueni-ovtchinnikov @KrisThielemans why create_uniform_copy is only available for PET? https://github.com/CCPPETMR/SIRF/blob/master/src/xSTIR/pSTIR/STIR.py#L347

KrisThielemans commented 5 years ago

we need it for an initialiser for many (traditional) PET algorithms. I guess it could be useful for MR as well, although I assume they will often start from a zero image. But I think good idea to move this up in the hierarchy. Could add a method in DataContainer that just implements it in terms of clone and then fill. Can still be overloaded for STIR.

(sounds a bit like we don't need it though, but I guess hard to get rid of it now)

paskino commented 5 years ago

We need it in CIL and it'd be useful to have it for MR as well. OK for the hint on how to implement it in DataContainer.

paskino commented 5 years ago

For STIR DataContainers there is no problem as the method create_uniform_copy exists. For Gadgetron

I don't quite know image templates exist
are we sure that the following (added to SIRF DataContainer) will work on templates? With PET it won't.
```
def get_uniform_copy(self, value=1.0):
y = self.clone()
y.fill(value)
return y
```

paskino commented 5 years ago

In CIL we have the class Operator which is equivalent to SIRF's AcquisitionModel. An Operator maps a space X -> Y. I.e. the image space to the acquisition space. If linear the Operator can map Y -> X.

We want to be able to query the Operator/AcquisitionModel about its space X and Y. For this I added 2 methods:

def range_geometry(self):
    return Y
def domain_geometry(self):
    return X

where I loosely call as Y and X the templates (or geometry in CIL). As of https://github.com/CCPPETMR/SIRF/pull/237/commits/a4cc4ec5e4e202730c2915584081216606508486 I save the AcquisitionData and ImageData used to set up the AcquisitionModel (PET and MRI) into instance variables. @evgueni-ovtchinnikov, there should be a way to retrieve the templates from the underlying C structures.

KrisThielemans commented 5 years ago

non-trivial. For ImageData we have the geometric info stuff but that is still incomplete as that doesn't capture time/repetition/etc. For AcquisitionData we don't have anything at all (aside from the complete data).

Not sure if we have the relevant info for Gadgetron.

paskino commented 5 years ago

closed by https://github.com/CCPPETMR/SIRF/commit/85c94c184c6b2d41387c217cedd2b9fe5913288a

SyneRBI / SIRF

get data container from acquisition model #215