Patlak \int_0_t Cp(t) wrongly computed in PatlakPlot.cxx

[AnderBiguri]

The following piece of code computes the Concentration of plasma in blood and its integral for the time frames used for the Patlak plot fit:

https://github.com/UCL/STIR/blob/7b4c27470ad175a80f29e9aabe945b6831a4f159/src/modelling_buildblock/PatlakPlot.cxx#L154-L167

It does so by integrating the IF over the time frame, and adds it to patlak_array[2][sample_num].

The integral part, it does by taking this same value, and multiplying by the frame duration and adding it to patlak_array[1][sample_num].

However, I think this last one is missing a big part. This being a Patlak plot, the first frame input will not be at t=0, but much later, lets call it t1. However, this code does not start integrating Cp(t) until t1. For an input where all time dynamic time frames are used (i.e. starting frame := 1 in the parameter file) and a starting frame that starts at t1, the first for loop in the code above simply does not run. sum_value=0 by the time the second loop starts, and thus the integral of Cp(t) is only ever computed from t1. But, as far as I know, we need the full integral of Cp(t), starting from zero.

I can verify numerically this is happening too, for a particular IF and a time frame length of 210, this is what the values of patlak_array are (after decay un-correction):

Cp(t)
1126.32 998.855 895.152 808.298 734.225 670.127 614.032 564.431 520.19 480.454 444.478 411.792 381.927 354.533 329.34 
int_Cp(t)
236527 441117 619456 775659 912891 1.03366e+06 1.14002e+06 1.23363e+06 1.3159e+06 1.38803e+06 1.45103e+06 1.50579e+06 1.55308e+06 1.59359e+06 1.62792e+06 

236527/210=1126.32

@KrisThielemans I'd appreciate a second pair of eyes to verify I'm not mad.

Should the integral not be from t=0? Or is this quite wrong?

[KrisThielemans]

Seems to be wrong indeed. However, in the other version of this function (i.e. PatlakPlot::get_model_matrix(const PlasmaData& plasma_data,const TimeFrameDefinitions& time_frame_definitions,const unsigned int starting_frame) ), it seems to be ok at first glance https://github.com/UCL/STIR/blob/7b4c27470ad175a80f29e9aabe945b6831a4f159/src/modelling_buildblock/PatlakPlot.cxx#L101-L104

I had noted already that these 2 functions are almost duplicates. I think only one of them is used. Let's hope it's the correct one... (then we can just remove the other)

[AnderBiguri]

Hum, I dont see any difference there. If the starting_frame is 1, then that code also does not run. I am quite sure we are using the one I linked as thats where I got the numerical values I posted.

[KrisThielemans]

I am quite sure we are using the one I linked as thats where I got the numerical values I posted.

ok. I see now that those lines were there anyway https://github.com/UCL/STIR/blob/19d349a12f133e84b346a5f5aaa99b6e043d77b9/src/modelling_buildblock/PatlakPlot.cxx#L157-L158

If the starting_frame is 1, then that code also does not run.

yes, but I guess in "normal" usage starting_frame would be something like 20. i.e. you have a PET study with a whole bunch of frames, but we'll only start Patlak for somewhere latish. Recall that in the code, plasma data is first put into same time frame definitions as images https://github.com/UCL/STIR/blob/c2d4cecd8bb5bd0b0e639790a0d307ad926fb0b3/src/modelling_buildblock/PatlakPlot.cxx#L409

This is why we need to create a bunch of dummy 0 images for the first 22 frames in the example. In fact, you'd really not only 1 dummy frame of the whole first 15min (or whatever you choose). I think it'd work then.

This is bad design of course. The whole integral thing should probably be rewritten.

[AnderBiguri]

Which I have not been doing, of course! ah!

In any case, I'll try to think a way to put the integral in there. We have the time and we have the data, there must be a relatively easy way of adding the integral without the need of an entire image.