STIR
STIR copied to clipboard
Add weights to Patlak linear regression
Weights can be turned on or off in par file, assume Poisson distribution
- This PR assumes data is NOT decay corrected.
@KrisThielemans Sure. Want to add some sort of check for decay correction? I don't remember well, but isn't one of the input parameters if the data is decay corrected?
interesting. Checking through PatlakPlot.cxx
there is no check for decay corrected images. All decay stuff is for the input-function.
Digging a bit deeper:
PatlakPlot.h
doxygen says that images are "in decaying counts", but it isn't checked. In fact it can't be checked as DynamicDiscretisedDensity
doesn't have a get_if_decay_corrected
(or simply decay_corrected()
).
Maybe you could also add that in? Should be a 5 min job? (hmmm...)
By the way, there's some outdated comments in the utility here. Once we remove that, we should also remove this stuff.
@KrisThielemans if I add to DynamicDiscretisedDensity
decay_corrected()
, this doesn't really solve much of a problem right? because we don't fill this nor use it anywhere. So I propose we just keep assuming everything is not decay corrected for now, adding the method to STIR that would check this and act accordingly seems like something that its its own PR, not just applicable to the Patlak stuff.
@AnderBiguri I've updated the doc. Have a look if it makes sense please.
Note that we have the capability with _in_total_cnt
to allow input in "activity" images (as opposed to usual STIR "counts"). I guess the Poisson weights would need modification then? (Maybe not because of the division with the cst of the matrix). could you check? If it gets too confusing, we could just error out for now...
@KrisThielemans when the input is in "activity", is it decay corrected too? My gut says that it would make sense for it to be called activity only if its decay corrected, but if it can be non-decay corrected, then its just a multiplication by a constant, right? If that is the case, even if it the weight was not taking the image values into account, it would be the same.
But indeed, the weights are reading the image values, so they are pixel_val/(\int Cp)^2
. I believe this equation only changes if there is decay correction, otherwise its just as is, irrelevant of pixel_val
units, right?
you're right that "activity" should be decay correct. We can have 3 cases:
- counts
- counts / frame_duration (not sure if there's an accepted name for this, but we could call it "non-decay-corrected activity" or similar)
- activity = counts / frame_duration * decay_factor
We currently say that we don't support decay corrected data, so I guess the last one is out. Only "counts" follow Poisson stats.
Regarding factors, global factors don't matter, as long as they're constant for every frame. However, that wouldn't be the case for any of these.
I think I wrote the eqs down somewhere but didn't get into the division by model-matrix. As this in itself has frame_duration and decay_factor contributions, this could be somewhat tricky, but then again for the weight calculation, it just seems a constant.
I have no time to check this better now. sorry. again, if it gets hairy, just throw an error and support only what you know is correct