AcquisitionModel::range_geometry() does not error if set_up with a template

[AnderBiguri]

If you do AcquisitionModel.set_up(AcquisitionData,ImageData) but AcquisitionData has been loaded from a template (I suppose this means there is no actual sinograms loaded), and then do

    x = acq_model.domain_geometry()
    y = acq_model.range_geometry()
    x.allocate(value = 'random')
    y_hat = acq_model.direct(x)
    y.allocate( value = 'random')
    x_hat = acq_model.adjoint(y)   # This fails 

The call to adjoint will fail (at least in NiftyPET). @rijobro says this may have to do with acq_model.range_geometry()

[KrisThielemans]

I'm fairly sure this doesn't have anything to do with the fact that set_up has been called or not. (set_up is fine for an AcquisitionModel is fine with the fact that the template is "empty" (i.e. only header and hence geometry, but pointing to an empty data file)).

y.allocate is trying to write to the object it was constructed with. I guess it makes sense, but it's still surprising if it points to an existing file.

maybe we should just document it?

[AnderBiguri]

Shouldn't functions check if the input parameters are valid? acq_model.adjoint(y) fails with something obscure as a message. if y is not valid, for whichever reason, should not that be checked somewhere?

[evgueni-ovtchinnikov]

y.allocate is trying to write to the object it was constructed with

No, it is not, it writes to a copy:

    def allocate(self, value=0, **kwargs):
        '''Alias to get_uniform_copy

        CIL/SIRF compatibility
        '''
        if value in ['random', 'random_int']:
            out = self.get_uniform_copy()
            shape = out.as_array().shape
            seed = kwargs.get('seed', None)
            if seed is not None:
                numpy.random.seed(seed) 
            if value == 'random':
                out.fill(numpy.random.random_sample(shape))
            elif value == 'random_int':
                max_value = kwargs.get('max_value', 100)
                out.fill(numpy.random.randint(max_value,size=shape))
        else:
            out = self.get_uniform_copy(value)
        return out

@AnderBiguri: So you should work with the returned object.

[KrisThielemans]

well, then the error comes from somewhere else. can you make a minimal example of the error. you don't need acq_model at all as far as I can see.

[AnderBiguri]

@KrisThielemans Following @evgueni-ovtchinnikov 's advice, if I do:

x = acq_model.domain_geometry()
y = acq_model.range_geometry()
x = x.allocate(value = 'random')
y_hat = acq_model.direct(x)
y = y.allocate( value = 'random')
x_hat = acq_model.adjoint(y)   # This does not fail

The error is gone and execution is correct. I assume I am way too late to make any impact on this, but should not allocate then have a different name?

[KrisThielemans]

now I'm even more confused. y.allocate returns a new object, and doesn't modify the existing?

allocate, direct etc names are fixed by CIL unfortunately.

[AnderBiguri]

While thanks to @evgueni-ovtchinnikov (thanks :D ) I can now run this code, I insist: if y is then not filled/allocated, should not adjoint catch that with a reasonable message?

[KrisThielemans]

I'll quote your error from https://github.com/UCL/STIR/pull/355

6: Traceback (most recent call last):
6:   File "/home/anderbiguri/.local/lib/python2.7/site-packages/nose/case.py", line 197, in runTest
6:     self.test(*self.arg)
6:   File "/home/anderbiguri/SIRF/build/sources/SIRF/src/xSTIR/pSTIR/tests/tests_niftypet.py", line 84, in test_main
6:     x_hat = acq_model.adjoint(y)
6:   File "/home/anderbiguri/SIRF/build/INSTALL/python/sirf/STIR.py", line 1261, in adjoint
6:     num_subsets = num_subsets)
6:   File "/home/anderbiguri/SIRF/build/INSTALL/python/sirf/STIR.py", line 1226, in backward
6:     check_status(image.handle)
6:   File "/home/anderbiguri/SIRF/build/INSTALL/python/sirf/Utilities.py", line 374, in check_status
6:     raise error(errorMsg)
6: error: ??? "'ProjDataFromStream: error reading data' exception caught at line 510 of /home/anderbiguri/SIRF/build/sources/SIRF/src/xSTIR/cSTIR/cstir.cpp; the reconstruction engine output may provide more information"

Python error reporting is always tedious, but if you look at the last one, it essentially says that STIR cannot read the data.

Error reporting in SIRF is very hard, as it doesn't necessarily know what the engine is going to do, but I cannot see SIRF testing at every function call if an object is actually readable or not. (In fact, SIRF cannot do that check anyway).

Feel free to suggest better wording for the STIR error and the SIRF error :-)

[AnderBiguri]

Well, my though when I got this error the first time was:

error:??? "'ProjDataFromStream: error reading data'

error: ????

The engine does not understand the error

ProjDataFromStream

I am using python, I have no idea what this object is

error reading data

Reading data from where? I have inputted the data.

So the engine gets an error that does not understand while doing something completely different from what I wanted to do.

I guess I am more used to the way MATLAB works assert(isempty(y), "Your data is empty"). I understand this is much more complex than that. I will just need to get used to this, as changing any of these messages would likely mean changing it everywhere.

[KrisThielemans]

I'm genuinely asking for suggestions on how to handle this. User experience is important. But it needs to be feasible.

I totally agree that the SIRF user will not know what a STIR object is. It probably wouldn't be impossible though to make this a little bit clearer. Original

error: ??? "'ProjDataFromStream: error reading data' exception caught at line 510 of /home/anderbiguri/SIRF/build/sources/SIRF/src/xSTIR/cSTIR/cstir.cpp; the reconstruction engine output may provide more information"

I don't think the line info is useful at all for the user(but it is for the developer I guess).

error: STIR threw the following exception: "'ProjDataFromStream: error reading data'

better? Feasible @evgueni-ovtchinnikov ?

I'm not so sure that it's unreasonable to have a "reading" error in acq_mode.adjoint(y). After all, it'll have to read from y. no?

[KrisThielemans]

STIR's message could potentially be improved (With some work) by including filenames, but that wouldn't help in most cases (as it's a ProjDataInMemory).

[AnderBiguri]

@KrisThielemans ah yes, sorry, I will reopen this then. This is much beyond my capabilities right now, as I can share my expectations, but for me STIR is a black box for now.

[paskino]

A brief resume on the CIL/SIRF integration:

Object	SIRF	CIL
volume	ImageData	ImageData
"sinogram"	AcquisitionData	AcquisitionData
volume metadata	ImageData	ImageGeometry
"sinogram" metadata	AcquisitionData	AcquisitionGeometry

In CIL we have 2 classes that contain the metadata for the DataContainer. Our specialised ImageData and AcquisitionData have inside a geometry object describing it. These geometries can allocate their respective ImageData or AcquisitionData.

In SIRF the metadata is held directly by the data itself and there isn't a geometry as such. A similar concept is a template where there is no data associated with the DataContainer.

The weird error you are getting is exactly due to this:

x = acq_model.domain_geometry()
x.allocate(value='random')

acq_model.direct(x)

1. x is a template, which doesn't contain data;
2. allocate returns a new object based on the template but with data, but it's discarded
3. calling direct on the template will fail with that strange error which tries to say that there's no data associated with the ImageData.

That error message should be made more clear indeed. However, the problem with the code above comes from 3 facts:

1. in SIRF there is no class to hold meta data of AcquisitionData or ImageData. These are held in the same classes.
2. The template does contain no data but it is of type AcquisitionData or ImageData giving the impression that it does hold the data
3. you (probably) find improper the use of a method allocate to allocate a different object

The only method @evgueni-ovtchinnikov and I found to allocate the data in a template is to call get_uniform_copy, which does return a new object rather than allocating the data in the template.

The code proposed by @evgueni-ovtchinnikov is the one that I would consider correct, except that I'd use

x = acq_model.domain_geometry().allocate(value='random')
y = acq_model.range_geometry().allocate(value='random')
y_hat = acq_model.direct(x)
x_hat = acq_model.adjoint(y)

Hopefully this clarifies a bit the matter.

[AnderBiguri]

@paskino much clearer indeed, thanks! definetly your way of using allocate makes it more clear what it does, I suppose I got confused due to the syntax I used.

[paskino]

BTW @KrisThielemans got it right at the first reply!!!

(set_up is fine for an AcquisitionModel is fine with the fact that the template is "empty" (i.e. only header and hence geometry, but pointing to an empty data file)).

However, he assumed that allocate would allocate the data in the template object rather than allocate a different object.

y.allocate is trying to write to the object it was constructed with. I guess it makes sense, but it's still surprising if it points to an existing file.

I reiterate, this behaviour is intended in CIL as the metadata is a different object with respect to the DataContainer, be it AcquisitionData or ImageData.

In SIRF, this may not be so clear as the metadata are held in an empty-of-data AcquisitionData or ImageData, which do not behave like DataContainers as they don't have the associated data. I find this unconventional and error prone as (I think) there's no way to know if the object is a template or not.

[evgueni-ovtchinnikov]

error: STIR threw the following exception:

ProjDataFromStream: error reading data

better? Feasible?

Feasible, but hardly better.

I don't think the line info is useful at all for the user (but it is for the developer I guess).

...but it is for the developer I guess: precisely.

The user is unlikely to fix the reported error. They have to send the error message to the developer, who would then know where to look first and might come back with possible reasons (reproducing the error situation may not be easy - user data and/or user OS may not be available for the developer etc.).

[AnderBiguri]

@paskino I wonder if there is something you can do inside allocate.

The code takes a copy: out = self.get_uniform_copy() Which will be empty if it is a template. However, then does: shape = out.as_array().shape out.fill(numpy.random.random_sample(shape))

So either fill will dump unused input values silently, or shape is [0,0,0,0] or something like that. Starting from the assumption that allocate should either return an allocated output or fail (to be consistent with the usage of allocate in programming in general), can we not just check this and error/warning (Warning/Error: Allocating an empty type(self).__name__)? Or alternatively, if the DataContainer is a template, it should have the information of shape somewhere else than in as_array().shape presumably (voxels ?). Can we not use that to allocate, thus if you allocate a template, it returns a full version of it?

@evgueni-ovtchinnikov from my limited experience with TIGRE this approach is quite dangerous. If the user needs the developer to understand things like wrong input data or wrong usage of a function, you either need to make sure all devs will be ready forever to answer all simple questions about the usage continuously (like now, with me hehe), or accept that people will stop using the software because they don't understand how to use it/what is wrong with the code. I assume the second is not what you want with SIRF, and the first one is so annoying that currently when I add things to TIGRE are generally to make sure people send me less emails :D . But again, I am the new one here, and you guys have been doing this amazing work for quite long, so I understand my input here is (and should be) limited. Just my two cents :)

[paskino]

@AnderBiguri get_uniform_copy uses just the metadata to allocate the data and fill it with a constant value, default being 1.0. Also it accepts the parameter value with which you can specify what value you want it filled with.

So if self is a template or DataContainer with data nothing changes and you'll get what you want.

WRT shape I'm more than happy not to use as_array().shape if there's a SIRF-esque way to do that without going through numpy.