2nd - nested Patlak and generalized Patlak OSMAPOSL (POSMAPOSL & GPOSMAPOSL) and nested Richardson-Lucy Motion Deconvolution Image Reconstruction (RLMCIR) by NikEfth · Pull Request #1736 · UCL/STIR

NikEfth · 2026-08-02T11:53:02Z

Changes in this pull request

This is the second attempt to introduce the generalized Patlak code and nested Richardson-Lucy Motion Deconvolution image reconstruction by @nkarakatsanis

Testing performed

Related issues

Checklist before requesting a review

[] I have performed a self-review of my code
[] I have added docstrings/doxygen in line with the guidance in the developer guide
[] I have implemented unit tests that cover any new or modified functionality (if applicable)
[] The code builds and runs on my machine
[] documentation/release_XXX.md has been updated with any functionality change (if applicable)

Contribution Notes

Please tick the following:

The content of this Pull Request (the Contribution) is intentionally submitted for inclusion in STIR (the Work) under the terms and conditions of the Apache-2.0 License.
I (or my institution) have signed the STIR Contribution License Agreement (not required for small changes).

codacy-production · 2026-08-02T11:54:19Z

Not up to standards ⛔

🔴 Issues 1 medium

Alerts:
⚠ 1 issue (≤ 0 issues of at least minor severity)

Results:
1 new issue

Category Results

BestPractice 1 medium

View in Codacy

🟢 Metrics 266 complexity · 128 duplication

Metric Results

Complexity 266

Duplication 128

View in Codacy

_{NEW Get contextual insights on your PRs based on Codacy's metrics, along with PR and Jira context, without leaving GitHub. Enable AI reviewer}
_{TIP This summary will be updated as you push new changes.}

NikEfth · 2026-08-10T04:55:36Z

I have noticed an important asymmetry with the parametric images IO. The problem is that KineticParameters is a template and that makes reading the image very hard as no class in the Input paradigms knows if the images will have 2 or 3 parameters beforehand.

@KrisThielemans I think we should make the kineticParameters non template in a new PR and merge here.

NikEfth · 2026-08-10T05:53:48Z

Needs #1744

KrisThielemans · 2026-08-10T07:21:59Z

I have noticed an important asymmetry with the parametric images IO. The problem is that KineticParameters is a template and that makes reading the image very hard as no class in the Input paradigms knows if the images will have 2 or 3 parameters beforehand.

@KrisThielemans I think we should make the kineticParameters non template in a new PR and merge here.

At the time, we made KineticParameters thinking it would speed-up things. Now I'm not so sure anymore. It causes dificulties, and it might be a good idea to remove it.

However,

not sure if it an easy solution for reading (you need to know how many parameters to read anyway)
if at all possible, I would not do that in this PR. It just means it's never going to be merged...

KrisThielemans

Some minor house-keeping comments for now, not a review.

copyright/headers seems to need changes in many places
no output to cout
seems you still need a precommit run --all. However, need to be careful with ownership then. Ideally, this is done as @nkarakatsanis before @NikEfth major changes.

KrisThielemans · 2026-08-10T07:12:54Z

@@ -0,0 +1,193 @@
+//
+/*
+    Copyright (C) 2006 - 2011, Hammersmith Imanet Ltd


surely incorrect

KrisThielemans · 2026-08-10T07:14:35Z

+  - The dynamic images will always be in decaying counts.
+  - The plasma data is assumed to be in decaying counts.
+
+  \todo Should be derived from LinearModels, but when non-linear models will be introduced, as well.


NikEfth · 2026-08-10T08:54:42Z

seems you still need a precommit run --all. However, need to be careful with ownership then. Ideally, this is done as @nkarakatsanis before @NikEfth major changes.

The first was ran as @nkarakatsanis

KrisThielemans · 2026-08-10T09:01:30Z

Great (run -all?). We can squash some of your commits (now or later) then.

NikEfth · 2026-08-10T09:16:57Z

not sure if it an easy solution for reading (you need to know how many parameters to read anyway)

if at all possible, I would not do that in this PR. It just means it's never going to be merged...

We cannot really ignore this because extract_single_images_from_parametric_image cannot read the 3-parametric images, and we cannot test easily in run_tests_modelling. Unless we add a different method...

KrisThielemans · 2026-08-10T10:26:09Z

I suggest to do something like this in ParametricDiscretisedDensity.h

  /! Convenience typedef for base-type of Cartesian Voxelised Parametric Images with just two parameters
  typedef VoxelsOnCartesianGrid<KineticParameters<2, float>> Parametric2VoxelsOnCartesianGridBaseType;
  //! Convenience typedef for Cartesian Voxelised Parametric Images with just two parameters
  typedef ParametricDiscretisedDensity<ParametricVoxels2OnCartesianGridBaseType> Parametric2VoxelsOnCartesianGrid;
  // same for 3
  ...
  //! type for 2 parameters, we need this for backwards compatibility
  using ParametricVoxelsOnCartesianGrid = Parametric2VoxelsOnCartesianGrid;

and in extract...

   auto param2_im_sptr(read_from_file<Parametric2VoxelsOnCartesianGrid>(argv[2]));
    if (param2_im_sptr)
       return call_function_that_does_the_extraction(param2_im_sptr);
   auto param3_im_sptr(read_from_file<Parametric2VoxelsOnCartesianGrid>(argv[2]));
    if (param3_im_sptr)
       return call_function_that_does_the_extraction(param2_im_sptr);
    error("Failed to read parametric image (" + std::string(argv[2]) + ").");

KrisThielemans · 2026-08-10T12:35:26Z

Not good. if you do a rebase -i, git forgets about merges. This will create major problems.

Either use git rebase --rebase-merges (which might work and keep original hashes, not sure), or rebase on top of master now such that you get necessray SWIG fixes (and stop merging with master). Please do this soon, as we'll get into a nightmare again otherwise.

BTW, please stop CI from running if we know it will fail.

NikEfth · 2026-08-10T17:08:40Z

git diff 1602403 3b86ead gives nothing. I think it is alright.

NikEfth · 2026-08-10T17:32:37Z

Not good. if you do a rebase -i, git forgets about merges. This will create major problems.
Either use git rebase --rebase-merges (which might work and keep original hashes, not sure), or rebase on top of master now such that you get necessray SWIG fixes

I think there not way around. The hashes will change.

(and stop merging with master).

:/ I am sorry for this old habits.

Please do this soon, as we'll get into a nightmare again otherwise.

Done.

BTW, please stop CI from running if we know it will fail.

KrisThielemans · 2026-08-11T06:56:45Z

Great. CI is now good!

@NikEfth I suggest to squash your current commits (maybe even all of them). There's little value in keeping them. I would even squash @nkarakatsanis 2 first commits (the actual code and the precommit): we can just keep the history for one. Any reason not to?

(cherry picked from commit d96f158) Run precommit.

1. Applied CODACY suggestions 2. Recon tests for NESTPOSMAPOSL * CODACY suggestions * NESTGPOSMAPOSL works. Test disabled for now because of inconstistaeny in parametric image reading precommit

codacy-production · 2026-08-12T18:57:14Z

Not up to standards ⛔

🔴 Issues 3 medium

Alerts:
⚠ 3 issues (≤ 0 issues of at least minor severity)

Results:
3 new issues

Category Results

BestPractice 3 medium

View in Codacy

🟢 Metrics 248 complexity · 103 duplication

Metric Results

Complexity 248

Duplication 103

View in Codacy

_{NEW Get contextual insights on your PRs based on Codacy's metrics, along with PR and Jira context, without leaving GitHub. Enable AI reviewer}
_{TIP This summary will be updated as you push new changes.}

forgotten copyright info fix copyright info

NikEfth · 2026-08-13T04:54:34Z

@KrisThielemans are you planning to squash the entire PR in the future? As it is now I am happy with the 3 commits

KrisThielemans · 2026-08-13T06:59:52Z

My plan is to merge. This will preserve authorship, which is better for everyone (and also for our use of git fame.

Current commits look great.

Of course, now I'll have to do a proper review...

NikEfth · 2026-08-13T08:07:36Z

The motion methods and IO and cerr stuff are not ready yet. I am sorry. A bit more

KrisThielemans · 2026-08-13T08:23:40Z

I'm not in a hurry to review!

…ad() checks if the number of parameters in the header matches the template. * Replaced Karakatsanis IO paths for 3-params images * read_interfile_parametric_image is templated on num_params * New typedefs Parametric3VoxelsOnCartesianGrid and Parametric2VoxelsOnCartesianGrid * write_basic_interfile is templated on num_params * PoissonNestedLogLikelihoodWithGeneralizedPatlakAndDynamicProjectionData and use ParseAndCreateFrom to create the Target. * cout and cerr replacement in many places

* set_matrix_in_total_frame_counts in PatlakPlot * Introduced NestedGeneralizedPatlakOSMAPOSL to initialize 3-params recon with a 2-params image. * Fixed the order of parameter names. * run_tests_modelling now works with NESTGPOSMAPOSL and compares the correct parameters * pre-commit

NikEfth · 2026-08-14T12:54:46Z

@KrisThielemans @nkarakatsanis
This is a good milestone.

[x] NESTPOSMAPOSL
[x] NESTGPOSMAPOSL

A lot of extra IO has been implemented for modernization and alignment with modern STIR.
Also now one can initialize NESTGPOSMAPOSL with a 2-params image.

I want to stress the lager number of "hacks" necessary to make this work with tempalated KineticParameters<num, float>. Not good.

NikEfth · 2026-08-14T16:08:26Z

+          double integral_step = cur_iter->get_plasma_counts_in_kBq()
+                                 * this->_plasma_frame_data.get_time_frame_definitions().get_duration(sample_num);
+
+          sum_value += 0.5 * integral_step;


@KrisThielemans @nkarakatsanis
If I understand this correctly. This moves the integral to the middle of the frame, rather than the beginning. This should break backwards compatibility. And more importantly test_modelling.cxx

Error : unequal values are 1996.12 and 1917.64. Check _model_array-1st column in ModelMatrix Error : unequal values are 2076.65 and 2005.63. Check _model_array-1st column in ModelMatrix Error : unequal values are 2141.73 and 2077.13. Check _model_array-1st column in ModelMatrix Error : unequal values are 2193.4 and 2134.47. Check _model_array-1st column in ModelMatrix Error : unequal values are 2233.58 and 2179.44. Check _model_array-1st column in ModelMatrix

(I prefer the new version)Any advice?

I'd have to review this more, but generally I don't trust KM code that assumes "mid-points" are ok. We measure counts in time frames. Strictly speaking, all integrals need to be computed over the time frame. However, I'm not sure if that's what the old code did (although I might have convinced Harry accordingly...).

But again, I don't know why @nkarakatsanis made that change, so let's see what he says!

Honestly I don't recall if I just followed an existing "KM convention" in STIR or simply made that choice. I have no strong feelings to change this to integration over the time frames. As long as we are explicit and consistent between plasma input function and TACs, I am fine with either approaches.

I think the GeneralizedPatlakPlot also follows this convention and is (cryptically) documented as:

std::cout << "\nCreating Generalized Model Matrix (Here printed in its transverse format)\n\n" << "NOTE1: It contains as many columns as the number of later frames participating in parameter estimation\n" << "It contains as many rows as the convolution points of the input function\n" << "+ 1 last row consisting of the plasma counts for the corresponding later frame\n" << "NOTE2: Last element of each column is the plasma counts for the corresponding later frame\n\n" << "The total number of frames are: " << this->_num_frames << "\n" << "The total number of complete plasma samples are: " << this->_complete_plasma_data.size() << "\n" << "The last frame middle time is : " << this->_last_frame_mid_time << "\n" << "The time shift in complete plasma samples is: " << this->_complete_plasma_data.get_time_shift() << "\n" << "The total number of convolution points + 1(one) more column are: " << this->_num_conv_params << "\n" << "First Column: plasma samples for frame 1 ... Last Column: plasma samples for last frame\n";

GeneralizedPatlakPlot has been designed to follow the same time frame convention with that implemented for PatlakPlot. Any changes will need to be implemented across both plots for consistent performance.

OK, I had to push an update although it is partial. This issue is resolved with the introduction of frame reference time := frame end | frame midpoint. We can discuss more

* PatlakPlot and GeneralizedPatlakPlot both derive from KineticModel, * Take the isotope half-life from the ExamInfo radionuclide * NESTGPOSMAPOSL is broken at the moment.

NikEfth · 2026-08-20T23:14:02Z

In order to avoid extenssive code duplication I did some major changes (ongoing).

PatlakPlot and GeneralizedPatlakPlot both derive from KineticModel
For initialization GeneralizedPatlakPlot holds a PatlakPlot member should the user want to use it. But can also initiliaze with an initial image from the par file.
frame reference time := frame end | frame midpoint allows the user to choose how to integrade the plasma values and is consisent between the two models.
The radionuclide info are taken from the ExamInfo. This is nof finished because the generate_image does not use ExamInfo, so this is deactivated for the moment.

The GeneralizedPatlakPlot has large portions commented out, for the moment. It needs further clean up, but this commit was getting a bit too big.

NikEfth added this to the v7.0 milestone Aug 2, 2026

NikEfth assigned NikEfth and unassigned NikEfth Aug 2, 2026

NikEfth marked this pull request as draft August 2, 2026 12:22

NikEfth mentioned this pull request Aug 2, 2026

nested Patlak and generalized Patlak OSMAPOSL (POSMAPOSL & GPOSMAPOSL) and nested Richardson-Lucy Motion Deconvolution Image Reconstruction (RLMCIR) #1685

Open

2 tasks

KrisThielemans requested changes Aug 10, 2026

View reviewed changes

NikEfth force-pushed the port/nkarakatsanis-current branch from 1602403 to 3b86ead Compare August 10, 2026 10:27

NikEfth force-pushed the port/nkarakatsanis-current branch from fb63adf to 09f50ea Compare August 10, 2026 17:30

nkarakatsanis and others added 2 commits August 12, 2026 21:54

Recover reconstruction and kinetic-modelling work

0d4858f

(cherry picked from commit d96f158) Run precommit.

Code modernisation and bringing up to speed with current master

e3c7d8f

1. Applied CODACY suggestions 2. Recon tests for NESTPOSMAPOSL * CODACY suggestions * NESTGPOSMAPOSL works. Test disabled for now because of inconstistaeny in parametric image reading precommit

NikEfth force-pushed the port/nkarakatsanis-current branch from 09f50ea to d6268e6 Compare August 12, 2026 18:56

Copyright messages

49ab80f

forgotten copyright info fix copyright info

NikEfth force-pushed the port/nkarakatsanis-current branch from d6268e6 to 49ab80f Compare August 12, 2026 18:58

NikEfth added 2 commits August 13, 2026 23:17

NikEfth marked this pull request as ready for review August 14, 2026 14:35

NikEfth commented Aug 14, 2026

View reviewed changes

NikEfth added 2 commits August 21, 2026 01:02

* Add configurable frame reference time to KineticModel

ae7592b

* PatlakPlot and GeneralizedPatlakPlot both derive from KineticModel, * Take the isotope half-life from the ExamInfo radionuclide * NESTGPOSMAPOSL is broken at the moment.

pre-commit

ae4bbb4

Uh oh!

Conversation

NikEfth commented Aug 2, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Changes in this pull request

Testing performed

Related issues

Checklist before requesting a review

Contribution Notes

Uh oh!

codacy-production Bot commented Aug 2, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Not up to standards ⛔

Uh oh!

NikEfth commented Aug 10, 2026

Uh oh!

NikEfth commented Aug 10, 2026

Uh oh!

KrisThielemans commented Aug 10, 2026

Uh oh!

KrisThielemans left a comment

Choose a reason for hiding this comment

Uh oh!

KrisThielemans Aug 10, 2026

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

KrisThielemans Aug 10, 2026

Choose a reason for hiding this comment

Uh oh!

NikEfth commented Aug 10, 2026

Uh oh!

KrisThielemans commented Aug 10, 2026

Uh oh!

NikEfth commented Aug 10, 2026

Uh oh!

KrisThielemans commented Aug 10, 2026

Uh oh!

KrisThielemans commented Aug 10, 2026

Uh oh!

NikEfth commented Aug 10, 2026

Uh oh!

NikEfth commented Aug 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

KrisThielemans commented Aug 11, 2026

Uh oh!

codacy-production Bot commented Aug 12, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Not up to standards ⛔

Uh oh!

NikEfth commented Aug 13, 2026

Uh oh!

KrisThielemans commented Aug 13, 2026

Uh oh!

NikEfth commented Aug 13, 2026

Uh oh!

KrisThielemans commented Aug 13, 2026

Uh oh!

NikEfth commented Aug 14, 2026

Uh oh!

NikEfth Aug 14, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

KrisThielemans Aug 14, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

nkarakatsanis Aug 14, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

NikEfth Aug 14, 2026

Choose a reason for hiding this comment

Uh oh!

nkarakatsanis Aug 14, 2026

Choose a reason for hiding this comment

NikEfth commented Aug 2, 2026 •

edited

Loading

codacy-production Bot commented Aug 2, 2026 •

edited

Loading

NikEfth commented Aug 10, 2026 •

edited

Loading

codacy-production Bot commented Aug 12, 2026 •

edited

Loading

NikEfth Aug 14, 2026 •

edited

Loading

KrisThielemans Aug 14, 2026 •

edited

Loading

nkarakatsanis Aug 14, 2026 •

edited

Loading

NikEfth commented Aug 20, 2026 •

edited

Loading