[Castor-users] Visualisation of a Line of the System Matrix

[M.R Teimoori]

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

On Thu, May 4, 2017 at 4:02 AM, STUTE Simon <simon.stute at cea.fr> wrote:

> 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.
>
>
>
> 1. 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.
>
>
>
> 2. 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:
>
>
>
> 1. 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).
>
>
>
> 2. 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
>
>
>
>
> ------------------------------
> *De :* M.R Teimoori [reza.teimoori.s at gmail.com]
> *Envoyé :* mercredi 3 mai 2017 01:57
> *À :* STUTE Simon; castor-users at lists.castor-project.org
> *Objet :* Re: [Castor-users] Visualisation of a Line of the System Matrix
>
> 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:
>
> - You probably see uniform images in iteration 1 and 2 and see your LOR(s)
> only after the 2nd iteration which puzzles me a bit,
> - Enabling multi-threading when only a few events are used directs the
> reconstruction code to an infinite loop.
>
> Thank you,
>
> Cheers,
> Reza
>
> On Tue, May 2, 2017 at 10:47 AM, Simon Stute <simon.stute at cea.fr> wrote:
>
>> 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
>>
>>
>> Le 28/04/2017 à 19:06, M.R Teimoori a écrit :
>>
>> 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
>>
>>
>>
>>
>> _______________________________________________
>> Castor-users mailing listCastor-users at lists.castor-project.orghttp://lists.castor-project.org/listinfo/castor-users
>>
>>
>>
>> _______________________________________________
>> Castor-users mailing list
>> Castor-users at lists.castor-project.org
>> http://lists.castor-project.org/listinfo/castor-users
>>
>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.castor-project.org/pipermail/castor-users/attachments/20170508/4229b2f0/attachment.html>