TOF sinogram-division leads to an error

[NicoleJurjew]

When I try to divide two TOF-sinograms, the operation itself works. When I try to do something with the object that is created, however, I get an error message. Here is what I did:

>>>acq_data = pet.AcquisitionData('true_emission.hs')
>>>print(acq_data.dimensions())
(33, 80, 399, 520)
>>>one_quot_enum = acq_data.get_uniform_copy(1)
>>>one_quot_denom = acq_data.get_uniform_copy(1)
>>>one_quot = one_quot_enum.divide(one_quot_denom)
>>>one_quot.as_array()

This leads to the error message attached.

[KrisThielemans]

Presumably this creates 3 temporary projdata (each with .hs and .s). Can you please check that headers are the same (aside from filename of course), and do ls -l *.s or similar to get file sizes

[NicoleJurjew]

The headers are all the same (example attached), but the file sizes differ. The two "get_uniform_copy(1)" sinograms are both 2.1 GB while the one that's created via division is only 1.1 GB big.

tmp_3_1648195191108.txt

[KrisThielemans]

ok. a bug.

I guess at the moment you will have to pass via as_array() to do the quotient.

[KrisThielemans]

@evgueni-ovtchinnikov any ideas? sizes are a bit weird (might even indicate it's filled with the wrong type as it's half the size)

[evgueni-ovtchinnikov]

@KrisThielemans the division is performed in PETAcquisitionData::binary_op using STIR's SegmentBySinogram iterators - do they work on TOF data with TOF dimension > 1?
@NicoleJurjew the current SIRF master cannot throw error at line 806 of cstir.cpp - you must have been using some older version of SIRF, please try to update and re-run your test

[KrisThielemans]

PETAcquisitionData::binary_op using STIR's SegmentBySinogram iterators - do they work on TOF data with TOF dimension > 1?

they should, in the sense that clone() and as_array() work normally, which presumably go via the same iterators.

Note that the test code above is with "in file" settings, not "in memory". Possibly this is something that @NicoleJurjew could test if the same problem occurs when using temp data in memory.

[evgueni-ovtchinnikov]

@NicoleJurjew can I get your true_emission files please?

[NicoleJurjew]

Thanks for your help, using "as_array" for the division I could solve my problem!

@evgueni-ovtchinnikov, @KrisThielemans , I sent a link to a OneDrive with the files via email. Please let me know if you can't access or need more information.

[evgueni-ovtchinnikov]

Quite a few operations on PETAcquisitionData are performed by loops like these:

float
PETAcquisitionData::norm() const
{
	double t = 0.0;
	for (int s = 0; s <= get_max_segment_num(); ++s)
	{
		SegmentBySinogram<float> seg = get_segment_by_sinogram(s);
		SegmentBySinogram<float>::full_iterator seg_iter;
		for (seg_iter = seg.begin_all(); seg_iter != seg.end_all();) {
			double r = *seg_iter++;
			t += r*r;
		}
		if (s != 0) {
			seg = get_segment_by_sinogram(-s);
			for (seg_iter = seg.begin_all(); seg_iter != seg.end_all();) {
				double r = *seg_iter++;
				t += r*r;
			}
		}
	}
	return sqrt((float)t);
}

As I discovered, these loops ignore the TOF dimension, which, in particular, gives reduced norm, and may be the reason for other trouble.

[KrisThielemans]

ok. this is a bit difficult.

The reason we have this is stir::ProjData doesn't have iterators at present as they'd be very inefficient in many cases. We only have them for ProjDataInMemory but that doesn't help (TOF sinograms can be too large to fit in memory anyway).

Easiest would be to have an extra loop over tof, but unfortunately it will need to be within an #ifdef STIR_TOF (or similar variable), as get_segment_by_sinogram doesn't have a TOF argument on STIR master of course.

A possible alternative is to make sure that all those loops are replaced by calls to equivalent STIR functionality. That is probably the most productive way of doing things, as they will be useful in STIR anyway. Of course, we will need to add the TOF loop there on the TOF branch, but that is quite natural.

@evgueni-ovtchinnikov which functionality are you missing in STIR with these loops?

[evgueni-ovtchinnikov]

@KrisThielemans the functionality we have problems at the moment is the data algebra: any algebraic operation on in-file PETAcquisitionData objects produces an object that cannot be cloned or as_array'ed - I get ProjDataFromStream: error reading data exception thrown. Surprisingly, the norm of such an object is computed relatively correctly (e.g. the norm of x + x is twice the norm of x, though both seem to be reduced to 1 TOF). It is also possible to use get_uniform_copy(), which creates a new object without pathologies (correct size and fully functional).

[KrisThielemans]

I had a quick look and see these loops in norm, dot, inv and binary_op (which seems only used in multiply, divide, maximum, minimum). These could be written in terms of a few basics like (using names that are like std` algorithms that work on sequences)

//! apply f (double to avoid  numerical precision) to accumulate
double ProjData::accumulate(double init, std::function<double(double, double)> f= std::plus<double>()) const;
//! apply f on each element of p and store in *this
void ProjData::transform(const ProjData p&, std::function<float(float)> f);
//! apply f elementwise on p1 and p2 and store in *this
void ProjData::transform(const ProjData& p1, const ProjData& p2, std::function<float(float,float)> f) ;

These wouldn't be too hard to write in STIR I believe. It'd allow cleaning-up some STIR code as well. (They'd have to be overloaded for ProjDataInMemory for efficiency.

Much better would be to write STIR iterators for ProjData that do the loop. That would need some thought, especially as overloading seems quite tricky.

What do you think?

[KrisThielemans]

@NicoleJurjew I think this can be closed since #1208?

[NicoleJurjew]

@KrisThielemans yes, I wasn't aware those were still floating around. Thanks!