Hi All,
I would like to visualize a sample line of the system matrix using CASToR. If I am not wrong, this should be feasible when you input a single coincidence data (corresponding to the desire line of the system matrix) and reconstruct an image using listmode EM with 1 iteration. However, CASToR returns a uniform image if the input has only the data for one coincidence. Why is this not possible? If anybody knows a workaround, please let me know.
Thank you for your help,
Cheers,
Reza
Hi again !
Could you send the geometry description of your scanner along with your tiny list-mode please ?
I will be able to test it quickly and find why you do not see nothing.
Cheers
Simon
Hello Simon,
Thank you for your replies to my two emails.
I’ve attached a file of a CASToR-formatted crystal lookup table for an arbitrary PET scanner called “aSamplePETScanner.lut” with its header “aSamplePETScanner.hscan”. Let’s assume that I’d like to visualize a sample line of the system matrix that corresponds to a coincidence between crystalIDs 6000 and 6142. The XYZ coordinates of the centroids of the crystals, as defined by CASToR’s coordinates system, are [-161.4 -11.0 -15.8] and [156.4 37.4 -12.6]. I’ve made a sample binary coincidence data file that has the information of this coincidence and set the time of the event to 1 ms. This file is called “aSampleLOR.Cdf” and is attached along with its header “aSampleLOR.Cdh”. No attenuation/scatter/random/normalization coefficient is defined in this overly simplified case.
I usually find visualizing a single row of the system matrix extremely helpful if I want to learn the effects of various reconstructions parameters on image quality. Today in my investigations, I learned that it actually is possible to see this coincidence only if no convolution is used. I’ve attached a 2D superimposed image of all image slices for the mentioned event to this email. Once I introduce a convolution (i.e. -conv … ), the images turn to be blank. So, I’m wondering how one can have the images of a single coincidence with some convolution parameters introduced.
Also, here are some side notes about what I observed about this discussion. They may be useful for other CASToR users too, if you use one or only a few coincidences:
Thank you,
Cheers,
Reza
aSampleLOR.Cdh (154 Bytes)
aSampleLOR.Cdf (12 Bytes)
aSamplePETScanner.lut (342 KB)
aSamplePETScanner.hscan (298 Bytes)
Hi Reza,
First, many thanks for experimenting with CASToR and pushing its boundaries!
I dug into the problems and found solutions.
When you are reconstructing from a list-mode, the whole sensitivity is taken into account. During an update (after a subset or an iteration of 1 subset in your case), if a voxel has never been visited while its sensitivity is not null, its update factor will be 0. In the standard case, we do not want that; in another subset (if more than 1 are used) the voxel may be updated so we do not want it to be stuck at a 0 value. So in the default configuration file of the MLEM algorithm (config/optimizer/MLEM.conf), we set a minimum update factor at 0.01. This explains why in your image you see your line only after a certain number of iterations. Indeed, as you do not have any calibration or so, the image value along your event should have a very low value. Then as the initial image value is 1, the update factor along your line is also limited to 0.01. So after a few iterations, the other voxels will have a sufficiently low value so as to let your event appear on the image.
And when you are using an image-based PSF, it spreads even more the line of response will lower values, explaining why you do not see it event after a few iterations.
To counter this, there are two easy solutions for you.
Change the parameters of the MLEM algorithm. Use the -help-opti option to get explanations about its 4 parameters. The first one is the image update, the second one is the threshold for the denominator value in the data-space update, and the third and fourth ones are the minimum and maximum allowed image update values. So if you use this option “-opti MLEM,1.,1.e-10,-1,-1”, there is no limitations about image updates. In this case all voxels no included in your single event will be 0 right after the first update, so your line will appear clearly, even with a PSF.
Use a single histogram event instead of a list-mode event. When reconstructing histogram data, it is assumed that all possible lines-of-response are included into the datafile. Then the sensitivity is computed on-the-fly with respect to this assumption. So with a single histogram event, the sensitivity image will only include your single event. Two consequences: your line will directly appear after one update (because all voxels with 0 sensitivity are set to 0); and there is no global sensitivity computation at the beginning of the program so this is faster.
In the help of the MLEM algorithm, it is said about the minimum image update parameter that a 0 or negative value will allow no minimum but still forbid a 0 update. This is not true and was part about an obsolete implementation prior to the first official release. So ignore this because 0 updates are allowed as soon as you parameterized the MLEM alogirhtm to do it.
About the infinite loop. Indeed, when you have more threads than the number of events, you get stuck into the loop over events. This is due to the implementation of the parallel reading of the events. Thank you for spotting that bug, we never tried this particular case! This is corrected for the next version.
Still, there are two workarounds for you now:
Use the option “-load 0” to specify that all events are read on-the-fly directly from the file, one by one (the events are not buffered into memory).
Restrict the number of threads used for projection computations to 1. You can still set a higher number of threads for image computations (convolutions, updates, etc). To do that, you can give two different number of threads to the thread option, as “-th 1,16” for example. This is explained in the computation help that can be obtained with option “-help-comp”. However, if you are using list-mode data, only one thread will be used for sensitivity computation so it will be slow, except if you provided the already computed sensitivity image using option “-sens”.
Thanks again for your feedback!
Cheers
Simon
Hello Simon,
I was away last week and only got a chance to see your email this morning. I just wanted to thank you and Thibaut for your clear answers.
Cheers,
Reza
