STIR icon indicating copy to clipboard operation
STIR copied to clipboard

Add weights to Patlak linear regression

Open AnderBiguri opened this issue 2 years ago • 6 comments

Weights can be turned on or off in par file, assume Poisson distribution

  • This PR assumes data is NOT decay corrected.

AnderBiguri avatar Sep 08 '21 13:09 AnderBiguri

@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?

AnderBiguri avatar Sep 15 '21 15:09 AnderBiguri

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 avatar Sep 15 '21 16:09 KrisThielemans

@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 avatar Sep 21 '21 16:09 AnderBiguri

@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 avatar Sep 22 '21 14:09 KrisThielemans

@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?

AnderBiguri avatar Sep 23 '21 10:09 AnderBiguri

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

KrisThielemans avatar Sep 23 '21 10:09 KrisThielemans