Enabling PET reconstruction from list mode data #1103

evgueni-ovtchinnikov · 2022-05-10T16:35:02Z

Fixes #134.

KrisThielemans · 2022-05-10T20:24:01Z

looks good to me. Now the objective function...

it should be virtually the same as the one for projection data. Main difference is that it only accepts a projection matrix (not a projector pair).

closes #1037

Edo lm recon

evgueni-ovtchinnikov · 2022-05-23T10:57:27Z

@NikEfth I increased my VM memory to 16384, and got past set up, but the job was killed (it appears, not long before it would finish):

sirfuser@vagrant:~$ cstir_test6
running test6.cpp...

WARNING: Interfile warning: 'default bin size (cm)' (0.208000) is expected to be 0.200000.


WARNING: Interfile warning: num_axial_blocks_per_bucket (35) is expected to be 7.


WARNING: Interfile warning: axial crystals per singles unit (8) is expected to be 56.


WARNING: Interfile warning: transaxial crystals per singles unit (8) is expected to be 32.


WARNING: Interfile warning: I have used all explicit settings for the scanner
	from the Interfile header, and remaining fields set from the
	UPENN_5rings model.


WARNING: Scanner UPENN_5rings: inconsistent transaxial block info

WARNING: Interfile parsing ended up with the following scanner:
Scanner parameters:= 
Scanner type := UPENN_5rings
Number of rings                          := 280
Number of detectors per ring             := 594
Inner ring diameter (cm)                 := 76.4
Average depth of interaction (cm)        := 0.7
Distance between rings (cm)              := 0.39655
Default bin size (cm)                    := 0.208
View offset (degrees)                    := 0
Maximum number of non-arc-corrected bins := 331
Default number of arc-corrected bins     := 331
Energy resolution := 0.109
Reference energy (in keV) := 511
Number of blocks per bucket in transaxial direction         := 4
Number of blocks per bucket in axial direction              := 35
Number of crystals per block in axial direction             := 8
Number of crystals per block in transaxial direction        := 8
Number of detector layers                                   := 1
Number of crystals per singles unit in axial direction      := 8
Number of crystals per singles unit in transaxial direction := 8
Scanner geometry (BlocksOnCylindrical/Cylindrical/Generic)  := Cylindrical
end scanner parameters:=


image dimensions: 559x200x200
Offset: (-208, -208, -0)
Spacing: (2.08, 2.08, 1.98275)
Size: (200, 200, 559)
Direction matrix: 
1, 0, 0
0, 1, 0
0, 0, -1

== ok

WARNING: WARNING: ProjMatrixByBinUsingRayTracing used for pixel size (in x,y) that is larger than the bin size.
Backprojecting with this matrix might have artefacts at views 0 and 90 degrees.


WARNING: WARNING: ProjMatrixByBinUsingRayTracing used for pixel size (in x,y) that is larger than the bin size.
Backprojecting with this matrix might have artefacts at views 0 and 90 degrees.

tangential LORs: 2
Setting cache path...
Cache path: /home/sirfuser/devel/buildVM/sources/SIRF/data/examples/TBPET/
PoissonLogLikelihoodWithLinearModelForMeanAndListModeData: Skipping input!
Setting scanner template...
Dummy LM file

WARNING: The default max ring difference is set based full Scanner geometry
Setting acquisition model...
Setting max ring diff. ...
Sensitivity images: /home/sirfuser/devel/buildVM/sources/SIRF/data/examples/TBPET/sens_%d.hv
Setting max cache size ...
Max cache: 1500000000
Setting objective function ...
Setting up the reconstruction, please wait ...

WARNING: WARNING: ProjMatrixByBinUsingRayTracing used for pixel size (in x,y) that is larger than the bin size.
Backprojecting with this matrix might have artefacts at views 0 and 90 degrees.


WARNING: WARNING: ProjMatrixByBinUsingRayTracing used for pixel size (in x,y) that is larger than the bin size.
Backprojecting with this matrix might have artefacts at views 0 and 90 degrees.


WARNING: WARNING: ProjMatrixByBinUsingRayTracing used for pixel size (in x,y) that is larger than the bin size.
Backprojecting with this matrix might have artefacts at views 0 and 90 degrees.


WARNING: We skip the check on balanced subsets and presume they are balanced!

WARNING: We skip the check on balanced subsets and presume they are balanced!
Reconstructor set up.
Starting loop with 10 threads
Starting loop with 10 threads
Starting loop with 10 threads
Starting loop with 10 threads
Starting loop with 10 threads
Starting loop with 10 threads
Starting loop with 10 threads
Starting loop with 10 threads
Starting loop with 10 threads
Killed
sirfuser@vagrant:~$

Would it be possible to shorten this test run time by quitting early?

NikEfth · 2022-05-23T11:03:22Z

Yes, can you try with this file?

…

On May 23, 2022, at 06:57, Evgueni Ovtchinnikov ***@***.***> wrote: @NikEfth <https://github.com/NikEfth> I increased my VM memory to 16384, and got past set up, but the job was killed (it appears, not long before it would finish): ***@***.***:~$ cstir_test6 running test6.cpp... WARNING: Interfile warning: 'default bin size (cm)' (0.208000) is expected to be 0.200000. WARNING: Interfile warning: num_axial_blocks_per_bucket (35) is expected to be 7. WARNING: Interfile warning: axial crystals per singles unit (8) is expected to be 56. WARNING: Interfile warning: transaxial crystals per singles unit (8) is expected to be 32. WARNING: Interfile warning: I have used all explicit settings for the scanner from the Interfile header, and remaining fields set from the UPENN_5rings model. WARNING: Scanner UPENN_5rings: inconsistent transaxial block info WARNING: Interfile parsing ended up with the following scanner: Scanner parameters:= Scanner type := UPENN_5rings Number of rings := 280 Number of detectors per ring := 594 Inner ring diameter (cm) := 76.4 Average depth of interaction (cm) := 0.7 Distance between rings (cm) := 0.39655 Default bin size (cm) := 0.208 View offset (degrees) := 0 Maximum number of non-arc-corrected bins := 331 Default number of arc-corrected bins := 331 Energy resolution := 0.109 Reference energy (in keV) := 511 Number of blocks per bucket in transaxial direction := 4 Number of blocks per bucket in axial direction := 35 Number of crystals per block in axial direction := 8 Number of crystals per block in transaxial direction := 8 Number of detector layers := 1 Number of crystals per singles unit in axial direction := 8 Number of crystals per singles unit in transaxial direction := 8 Scanner geometry (BlocksOnCylindrical/Cylindrical/Generic) := Cylindrical end scanner parameters:= image dimensions: 559x200x200 Offset: (-208, -208, -0) Spacing: (2.08, 2.08, 1.98275) Size: (200, 200, 559) Direction matrix: 1, 0, 0 0, 1, 0 0, 0, -1 == ok WARNING: WARNING: ProjMatrixByBinUsingRayTracing used for pixel size (in x,y) that is larger than the bin size. Backprojecting with this matrix might have artefacts at views 0 and 90 degrees. WARNING: WARNING: ProjMatrixByBinUsingRayTracing used for pixel size (in x,y) that is larger than the bin size. Backprojecting with this matrix might have artefacts at views 0 and 90 degrees. tangential LORs: 2 Setting cache path... Cache path: /home/sirfuser/devel/buildVM/sources/SIRF/data/examples/TBPET/ PoissonLogLikelihoodWithLinearModelForMeanAndListModeData: Skipping input! Setting scanner template... Dummy LM file WARNING: The default max ring difference is set based full Scanner geometry Setting acquisition model... Setting max ring diff. ... Sensitivity images: /home/sirfuser/devel/buildVM/sources/SIRF/data/examples/TBPET/sens_%d.hv Setting max cache size ... Max cache: 1500000000 Setting objective function ... Setting up the reconstruction, please wait ... WARNING: WARNING: ProjMatrixByBinUsingRayTracing used for pixel size (in x,y) that is larger than the bin size. Backprojecting with this matrix might have artefacts at views 0 and 90 degrees. WARNING: WARNING: ProjMatrixByBinUsingRayTracing used for pixel size (in x,y) that is larger than the bin size. Backprojecting with this matrix might have artefacts at views 0 and 90 degrees. WARNING: WARNING: ProjMatrixByBinUsingRayTracing used for pixel size (in x,y) that is larger than the bin size. Backprojecting with this matrix might have artefacts at views 0 and 90 degrees. WARNING: We skip the check on balanced subsets and presume they are balanced! WARNING: We skip the check on balanced subsets and presume they are balanced! Reconstructor set up. Starting loop with 10 threads Starting loop with 10 threads Starting loop with 10 threads Starting loop with 10 threads Starting loop with 10 threads Starting loop with 10 threads Starting loop with 10 threads Starting loop with 10 threads Starting loop with 10 threads Killed ***@***.***:~$ Would it be possible to shorten this test run time by quitting early? — Reply to this email directly, view it on GitHub <#1103 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/ACEUB7U2TXCVNVKYSQ4TD43VLNQCLANCNFSM5VSFWYIA>. You are receiving this because you were mentioned.

evgueni-ovtchinnikov · 2022-05-23T11:31:42Z

@NikEfth Thanks a lot, it works much faster now and does not get killed!

NikEfth · 2022-05-23T11:34:15Z

Sure, but the images you get should not make sense

evgueni-ovtchinnikov · 2022-05-23T11:37:00Z

@NikEfth not important at this stage

… data

KrisThielemans

some minor comments which should be safe to implement as far as I can see.

src/xSTIR/cSTIR/include/sirf/STIR/stir_data_containers.h

src/xSTIR/cSTIR/tests/run_test6.cpp

src/xSTIR/cSTIR/tests/test6.cpp

src/xSTIR/pSTIR/STIR.py

KrisThielemans · 2022-05-25T06:41:54Z

src/xSTIR/pSTIR/STIR.py

-    if likelihood_type == 'LinearModelForMean':
+    #if likelihood_type == 'LinearModelForMean':
+    if likelihood_type is None:
        obj_fun = PoissonLogLikelihoodWithLinearModelForMeanAndProjData()
-        obj_fun.set_acquisition_data(acq_data)
+        if acq_data is not None:
+            obj_fun.set_acquisition_data(acq_data)


not sure why we're changing this code here. Do we really need to? What happens if acq_data is None? It'd presumably crash. (Only the listmode case can handle this)

@KrisThielemans so far as I can see, Nikos' list mode objective function does not need acquisition data, so acq_data must be optional. We need to discuss this with Nikos.

Agreed. @NikEfth's listmode function can either use listmode data or the corresponding cache.

However, these lines are for the "old" case with projection data, which do need it.

ok, will make this function throw in the non-listmode case

presumably the old code did already, but if not, best to put it in, agreed.

I still think we can revert the default to what it as, such that we have minimal code change, but as you prefer.

KrisThielemans · 2024-01-18T15:54:13Z

However, then I don't know what to do for the reconstruction.

@evgueni-ovtchinnikov @paskino Introducing an almost empty class on top of DataContainer indeed allows me to use C++ RTTI, so this problem is resolved like this

SIRF/src/xSTIR/cSTIR/cstir_p.cpp

Lines 1056 to 1064 in 61e692d

    
           SPTR_FROM_HANDLE(ContainerBase, sptr_cont, hv); 
        
           if (auto sptr_ad = std::dynamic_pointer_cast<STIRAcquisitionData>(sptr_cont)) { 
        
           	recon.set_input_data(sptr_ad->data_sptr()); 
        
           } 
        
           else if (auto sptr_ld = std::dynamic_pointer_cast<ListmodeData>(sptr_cont)) { 
        
           	recon.set_input_data(sptr_ld->data_sptr()); 
        
           } 
        
           else 
        
           THROW("input_data needs to be either ListmodeData or AcquisitionData");

So I think we can go ahead with this approach.

I haven't checked the Python side yet. We did introduce a ScanData class in Python, but have now no correspondence in C++ anymore. But maybe we don't need it.

- derive ListmodeData from ContainerBase - add ListmodeData to C wrappers - add get_info() to STIRImageData and ListmodeData - expand STIRAcquisitionData::get_info() to add exam-info string - minimal documentation This now works fine through Python, but demos still need work.

KrisThielemans · 2024-02-01T09:52:13Z

things to do:

Python test
remove osem_lm_reconstruction.py example
update listmode_reconstruction.py example to use ListmodeData object as opposed to filename
change name of objective function

- add ListmodeData::acquisition_data_template() - let ListmodeToSinograms::set_input() cope with either filename or ListmodeData object

KrisThielemans · 2024-03-28T03:13:19Z

listmode_reconstruction.py example now works (i.e. runs through and produces an image. Not 100% sure if it is correct. Could be that the acquisition model isn't passed through completely.

ListmodeToSinograms.set_input() now accepts an ListmodeData object as well.

KrisThielemans · 2024-04-04T16:48:30Z

Best to remind all use_gpu lines from the demo, as it isn't supported.

KrisThielemans

Actually, the whole osem_lm_reconstruction.py needs to be removed. I think everything in it is in listmode_reconstruction.py, which is more complete. Some things on the osem demo are currently not supported anyway.

KrisThielemans · 2024-04-04T16:41:54Z

examples/Python/PET/reconstruct_from_listmode.py

@@ -23,7 +23,7 @@
  -C <cnts>, --counts=<cnts>   account for delay between injection and acquisition start by shifting interval to start when counts exceed given threshold.
  --visualisations             show visualisations
  --nifti                      save output as nifti
-  --gpu                        use gpu
+  --gpu                        use gpu (actually means: use NiftyPET projectors)


Please keep removed. The STIR listmode recon only supports projmatrixbybin at the moment.

do you mean go back to just use gpu?

No. I mean remove all occurrences of use_gpu and the relevant lines to use niftypet, as it won't work anyway. See also main comment

KrisThielemans · 2024-04-04T16:58:33Z

I think that later we should move the preparation of randoms, scatter etc to SIRF itself, i.e. taking most of what you wrote for the challenge. However, let's do that in a separate PR, and merge this one asap.

evgueni-ovtchinnikov · 2024-04-05T09:03:04Z

Still have trouble with listmode_reconstruction.py:

sirfuser@vagrant:~/devel/buildVM/sources/SIRF/examples/Python/PET$ python3 listmode_reconstruction.py 
Finding files in /home/sirfuser/devel/install/share/SIRF-3.6/data/examples/PET/mMR

INFO: CListModeDataECAT8_32bit: opening file /home/sirfuser/devel/install/share/SIRF-3.6/data/examples/PET/mMR/list.l

Processing time frame 1
3500000 events stored
Number of prompts stored in this time period : 3612640
Number of delayeds stored in this time period: 0
Last stored event was recorded before time-tick at 50 secs
Total number of counts (either prompts/trues/delayeds) stored: 3612640

This took 2.9s CPU time.
estimate_randoms: trying GEHDF5...
HDF5-DIAG: Error detected in HDF5 (1.10.7) thread 1:
  #000: ../../../src/H5F.c line 292 in H5Fis_hdf5(): unable to determine if file is accessible as HDF5
    major: File accessibility
    minor: Not an HDF5 file
  #001: ../../../src/H5Fint.c line 929 in H5F__is_hdf5(): unable to open file
    major: File accessibility
    minor: Unable to initialize object
  #002: ../../../src/H5FD.c line 741 in H5FD_open(): open failed
    major: Virtual File Layer
    minor: Unable to initialize object
  #003: ../../../src/H5FDsec2.c line 360 in H5FD__sec2_open(): unable to open file: name = 'UNKNOWN', errno = 2, error message = 'No such file or directory', flags = 0, o_flags = 0
    major: File accessibility
    minor: Unable to open file
not a GE HDF5 file. Using ML estimate from delayeds
processed frame 1
Last stored event was recorded after time-tick at 49.999 secs
Total number of prompts/trues/delayed stored: 410099

This took 0.4s CPU time.
Starting iteration 1	KL 247.329
Starting iteration 2	KL 18.2321
Starting iteration 3	KL 0.663053
Starting iteration 4	KL 0.0278851
Starting iteration 5	KL 0.00129212
Starting iteration 6	KL 4.15109e-05
Starting iteration 7	KL 2.26989e-06
Starting iteration 8	KL 3.85613e-08
Starting iteration 9	KL 4.10926e-09
Starting iteration 10	KL 5.79696e-10
CPU time 30.69 secs

WARNING: DiscretisedDensity does not contain any time frames. This might cause an error.
Reading manufacturer PET normalisation file from /home/sirfuser/devel/install/share/SIRF-3.6/data/examples/PET/mMR/norm.n.hdr
trying GEHDF5...
not a GE HDF5 file
applying attenuation (please wait, may take a while)...
estimating scatter (this will take a while!)
Killed

(see also our exchange with Nikos on 23 May 2022 above - back then reducing the size of data helped).

KrisThielemans · 2024-04-05T14:16:48Z

I guess this is a problem on your VM. You should probably give it more memory.

It's not clear where this happens. If it's during scatter estimation (as the above implies), reducing number of listmode events is not going to change anything, as the scatter estimate is a (full resolution) projdata.

In any case, the number of events used here is very small.

KrisThielemans · 2024-04-05T14:17:10Z

I think that later we should move the preparation of randoms, scatter etc to SIRF itself, i.e. taking most of what you wrote for the challenge. However, let's do that in a separate PR, and merge this one asap.

evgueni-ovtchinnikov · 2024-04-05T17:04:28Z

I guess this is a problem on your VM. You should probably give it more memory.

True. I have increased VM memory from 6GB to 12GB, and listmode_reconstruction.py now runs ok (takes a lot of time though).

Permission to merge?

KrisThielemans · 2024-04-05T17:14:06Z

Great. Did you remove the osem_lm_reconstruction.py? Otherwise good to go.

…econstruction.py

KrisThielemans · 2024-04-05T20:34:32Z

Oh. Don't forget to add something to the doc and changes.md

implemented basic-functionality listmode data class in C++ and Python

51fed6c

evgueni-ovtchinnikov and others added 14 commits May 11, 2022 19:34

added objective function type for lismode reconstruction

c4be9dd

introduced abstract base class PETScanData for all raw data objects

0063f78

[ci skip] interfaced PETScanData into Python

b143cfc

Minimal backbone for LM recons

b209f33

New test for LM data (test6)

31ef16c

fix

b38b35b

Merge branch 'HEAD' into lm-recon

036c7c2

Working test6

65f9fe0

Generalize the paths in test6

97f76cd

Minor renaming

26d36fc

simplify link library to Boost

db4479d

closes #1037

Merge branch 'master' into HEAD

debcd5b

remove dangling endif

b59524c

Merge pull request #1 from paskino/edo_lm-recon

f65e1de

Edo lm recon

evgueni-ovtchinnikov added 5 commits May 23, 2022 12:13

[ci skip] moved data folder TBPET out of SIRF repo

e4632db

removed unused stuff from cstir_test6.cpp

ea18b8d

restored backward compatibility of run_test6.cpp

7a33289

implemented C interface for list mode reconstruction objective function

7e0e299

started interfacing reconstruction from listmode into Python

47eb909

KrisThielemans added this to the v3.3 milestone May 24, 2022

implemented minimal Python interface for reconstruction from listmode…

08f0425

… data

KrisThielemans requested changes May 25, 2022

View reviewed changes

resolved some issues raised by Kris

4dd2de7

evgueni-ovtchinnikov and others added 5 commits January 19, 2024 13:05

attended to Codacy issues

27eec3d

attended to further Codacy issues

81f1389

renamed ListmodeData to STIRListmodeData for consistency

b98c6d0

add set_acquisition_model to listmode obj-fun

713ad10

KrisThielemans added 4 commits March 28, 2024 02:57

minor doc corrections

835b499

listmode recon improvements

9ce3b51

- add ListmodeData::acquisition_data_template() - let ListmodeToSinograms::set_input() cope with either filename or ListmodeData object

Merge remote-tracking branch 'origin/master' into lm-recon

f009fbf

minor fix in demo

bb42564

KrisThielemans and others added 3 commits March 28, 2024 14:41

fix storing filename of listmode file

56fbace

add SIRF_DATA_PATH to STIR C++ test

b384409

fixed bugs and issues found in lmrecontest2

95b3a24

KrisThielemans requested changes Apr 4, 2024

View reviewed changes

[ci skip] got rid of gpu stuff in reconstruct_from_listmode.py

bf361ff

[ci skip] removed osem_lm_reconstruction.py, superseded by listmode_r…

752c077

…econstruction.py

[ci skip] updated CHANGES.md

e4ac001

evgueni-ovtchinnikov merged commit 2c66faf into master Apr 8, 2024
1 check passed

evgueni-ovtchinnikov deleted the lm-recon branch April 8, 2024 11:11

KrisThielemans mentioned this pull request Apr 9, 2024

Lmrecontest2 #1240

Closed

6 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Enabling PET reconstruction from list mode data #1103

Enabling PET reconstruction from list mode data #1103

evgueni-ovtchinnikov commented May 10, 2022

KrisThielemans commented May 10, 2022

evgueni-ovtchinnikov commented May 23, 2022

NikEfth commented May 23, 2022 via email

evgueni-ovtchinnikov commented May 23, 2022

NikEfth commented May 23, 2022

evgueni-ovtchinnikov commented May 23, 2022

KrisThielemans left a comment

KrisThielemans May 25, 2022

evgueni-ovtchinnikov May 25, 2022

KrisThielemans May 25, 2022

evgueni-ovtchinnikov May 25, 2022

KrisThielemans May 25, 2022

KrisThielemans commented Jan 18, 2024

KrisThielemans commented Feb 1, 2024 •

edited

Loading

KrisThielemans commented Mar 28, 2024

KrisThielemans commented Apr 4, 2024

KrisThielemans left a comment

KrisThielemans Apr 4, 2024

evgueni-ovtchinnikov Apr 5, 2024

KrisThielemans Apr 5, 2024

KrisThielemans commented Apr 4, 2024

evgueni-ovtchinnikov commented Apr 5, 2024

KrisThielemans commented Apr 5, 2024

KrisThielemans commented Apr 5, 2024

evgueni-ovtchinnikov commented Apr 5, 2024

KrisThielemans commented Apr 5, 2024

KrisThielemans commented Apr 5, 2024

Enabling PET reconstruction from list mode data #1103

Enabling PET reconstruction from list mode data #1103

Conversation

evgueni-ovtchinnikov commented May 10, 2022

KrisThielemans commented May 10, 2022

evgueni-ovtchinnikov commented May 23, 2022

NikEfth commented May 23, 2022 via email

evgueni-ovtchinnikov commented May 23, 2022

NikEfth commented May 23, 2022

evgueni-ovtchinnikov commented May 23, 2022

KrisThielemans left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

KrisThielemans commented Jan 18, 2024

KrisThielemans commented Feb 1, 2024 • edited Loading

KrisThielemans commented Mar 28, 2024

KrisThielemans commented Apr 4, 2024

KrisThielemans left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

KrisThielemans commented Apr 4, 2024

evgueni-ovtchinnikov commented Apr 5, 2024

KrisThielemans commented Apr 5, 2024

KrisThielemans commented Apr 5, 2024

evgueni-ovtchinnikov commented Apr 5, 2024

KrisThielemans commented Apr 5, 2024

KrisThielemans commented Apr 5, 2024

KrisThielemans commented Feb 1, 2024 •

edited

Loading