List mode data and LUT

[ASK]

Dear Castor user, I am new to Castor and imaging in general. I am trying to produce an image from data collected from a PET scanner we made. I have list mode data, with channel IDs and TOF and a LUT for our scanner. I have taken only around 100K coincidence from our list mode data as the total coincidence is around 33 million and I figured less coincidence would result in faster processing speed and I just want to test if I could receive any decent images. I have currently, modified a preexisting geom scanner from the config folder to make a rough geometry for our scanner since I could not figure out how to input our LUT to make the scanner. I then used castor-datafileConversion with the 100K coincidence list-mode data and the roughly made scanner and created the cdf and cdh file I want to use TOF reconstruction but wasn’t sure how to activate it. I have tried many times to look over the general document specifically the section of CASToR input datafile format and saw that there is a TOF information flag but i am not sure how to activate it. I have tried reading the C++ code but could see anything in the datafileConversion for TOF. I also tried using vim to modify the cdh file and add in the TOF flag and other information. I then used castor-recon with 10 iteration and 16 subsets. The command runs and I obtain iteration images and a sensitivity image, but when I use amide to look at the images I see that only iteration 1 (the whiter/more dense white space image) and the sensitivity image have an image and the other iterations are blank/black image. My question is:

1. Why am I not seeing any images past iteration 1?I will upload the iteration 1 image and iteration 2 image(which is the same for iterations 2-10) and the sensitivity image.
2. How can I use my LUT to make a custom-made scanner and is there any way to verify if a scanner geometry is correct or not?

[tmerlin]

Hi,
Firstly I advice you to ignore TOF until you get everything right regarding the conversion of the detector IDs and the data corrections. If you want to set it up anyway, the flag to put in the datafile header are TOF information flag: (set to 1) and TOF resolution (ps): with the proper time resolution. The PET list-mode benchmark contains TOF data with quantized measurements if you need an example.

Regarding the amount of coincidences I would rather suggest to take at least several millions of coincidence, or even the whole data as the reconstruction time should not be that long with a regular computer and multi-threading enabled (actually for list-mode, the computation of sensitivity image could take most of the reconstruction time depending on the number of data channels in your system).

I believe what you see in the iteration 1 is the initial value of the image (0.01 for the voxels in the FOV by default with MLEM). If you don’t see anything in the following images, either this is a contrast issue or there is indeed no data which means probably something is wrong in the geometry or data conversion.
Looking at your sensitivity image, I think the geometry seems right and I would double check whether the data channels are associated to the correct detectors IDs.

There are some toolkits you could use to check the consistency of your data:

• castor-scannerLUTExplorer: To check the actual position of the IDs from the geom or LUT file. So this can be used to check if the IDs are located where you expect them to be. You can use it to print all the IDs positions or check if element by element. If you want to check your LUT, you need to give its header to the -sf option.

$ castor-scannerLUTExplorer 

Usage:  castor-scannerLUTExplorer  -df datafile.cdh  [settings]

This program can be used to explore a datafile and get some info about it. By default, it simply prints general information recovered
from the reader. If the '-e' option is supply, an element-by-element exploration will be performed: information about each element will be
displayed.
For CT scanners or SPECT cameras described by a .geom file, a datafile must also be provided in order to get information regarding the projection angles.

[Mandatory parameters]:
  -sf   scan_file      : Give the path to the header of a single lut file, or to a geom file.

[Options]:
  -e                   : Flag for scanner element by element exploration (will list the content of each element).
  -g                   : Flag for global exploration (geometric characteristics of all elements will be displayed).
  -o   path_out_file   : Path to an output file in which the data will be printed.
  -df  data_file       : For CT or SPECT system, using a .geom file as input. Give the path to the header of a datafile in order to get projection informations.
  -vb value            : Give the verbosity level, from 0 (no verbose) to 2 (default: 1)
  --help,-h,-help      : Print out this help page.

• castor-datafileExplorer: To check event by event the content of the datafile (associated IDs and correction data).

$ castor-datafileExplorer 

Usage:  castor-datafileExplorer  -df datafile.cdh  [settings]

This program can be used to explore a datafile and get some info about it. By default, it simply prints general information recovered
from the reader. If the '-e' option is supply, an event-by-event exploration will be performed: information about each event will be
displayed.

[Mandatory parameters]:
  -df   datafile.cdh   : Give the path to a single datafile.

[Options]:
  -e                   : Flag for event by event exploration (will list the content of each event).
  -i                   : Flag for interactive one by one event exploration when -e option is supplied.
  -g                   : Flag for global exploration (will not list the content of each event but will supply a summary).
  -vb value            : Give the verbosity level, from 0 (no verbose) to 2 (default: 1)
  --help,-h,-help      : Print out this help page

Hope this helps!

[ASK]

Thank you for the reply! Could you guide me on what parts of the datafileConversion.cc I would need to write. For the moment I have written

uint32_t nEvents =14000000;

// ===> Acquisition duration and start time in seconds
// ===> This should be recovered from the dataset metadata)
FLTNB start_time_sec = 0.;
FLTNB stop_time_sec = 100.;
FLTNB duration_sec = stop_time_sec - start_time_sec;

// ===> Scanner variables (blocs, modules, submodules, crystals etc…)
// ===> depending on the system.
uint32_t nRsectorsPerRing = 38;
uint32_t nb_crystals_trs = 8,
nb_crystals_axl = 8;

What else would I need to write so that my events contains the IDs related to my data?

[ASK]

I have tired erasing the
// ROOT objects/variables
int32_t rSectorID1 = 0, rSectorID2 = 0;
int32_t moduleID1 = 0, moduleID2 = 0;
int32_t crystalID1 = 0, crystalID2 = 0;

#ifdef CASTOR_ROOT
TApplication Tapp = new TApplication(“tapp”,0,0);
TTree Coincidences = new TTree;
TFile* Tfile_root = new TFile(path_to_input_file.c_str(),“READ”,“ROOT file with histograms”);
Coincidences = (TTree*)Tfile_root->Get(“Coincidences”);
nEvents += Coincidences->GetEntries();

//------------- ROOT branches
double_t time1=0; 
Coincidences->SetBranchAddress("time1",&time1);
Coincidences->SetBranchAddress("rsectorID1",&rSectorID1);
Coincidences->SetBranchAddress("moduleID1",&moduleID1);
Coincidences->SetBranchAddress("crystalID1",&crystalID1);

double_t time2=0;
Coincidences->SetBranchAddress("time2",&time2);
Coincidences->SetBranchAddress("rsectorID2",&rSectorID2);
Coincidences->SetBranchAddress("moduleID2",&moduleID2);
Coincidences->SetBranchAddress("crystalID2",&crystalID2); 

// Recover acquisition start time from the first event
Coincidences->GetEntry(0);

#endif
and
#ifdef CASTOR_ROOT
Coincidences->GetEntry(i);
time_ms = time1;
#endif

// ===> Compute CASToR indexation depending on the system layout
// ===> This directly depends on the type of scanner file used to
// ===> described the system geometry.
// ===> See documentation for more information
uint32_t crystals_trs_id1 = (uint32_t)(crystalID1/nb_crystals_trs);

uint32_t ringID1 = moduleID1*nb_crystals_axl 
                 + crystals_trs_id1;
                 
uint32_t crystals_trs_id2 = (uint32_t)(crystalID2/nb_crystals_trs);

uint32_t ringID2 = moduleID2*nb_crystals_axl 
                 + crystals_trs_id2;

// Loop on the LORs contained in the event
for(int l=0 ; l<nb_lines_in_event ; l++)
{
  // Recover the CASToR ID position of the two crystals according to
  // the LUT corresponding to this scanner
  castor_id1[l] = ringID1*nb_crystals_per_ring 
                + rSectorID1*nb_crystals_trs 
                + crystalID1%nb_crystals_trs ;
                
  castor_id2[l] = ringID2*nb_crystals_per_ring 
                + rSectorID2*nb_crystals_trs
                + crystalID2%nb_crystals_trs ;
}

since in the datafileConversion it says those lines “Must be erased when adjusting this code to the conversion of any dataset” but when I do this I get this error when running the castor-recon |***** iScannerPET::GetPositionsAndOrientations() → Crystal index 1 (17236080) out of range [0:9727] !

***** vProjector::Project() → A problem occurred while getting positions and orientations from scanner !

***** oProjectorManager::ComputeProjectionLine() → A problem occurred while forward projecting an event !

***** iIterativeAlgorithm::StepInnerLoopInsideSubsetLoop() → A problem occurred while computing the projection line !

-------------------------------------------------|

***** iIterativeAlgorithm::StepInnerLoopInsideSubsetLoop() → A problem occurred inside the parallel loop over events !

***** vAlgorithm::IterateCPU() → A problem occurred while calling StepInnerLoopInsideSubsetLoop() function !

***** castor-recon() → Error while performing the reconstruction

furthermore, using castor-datafileExplorer for every event ID it showing
iEventListPET::Describe() → Display contents

Time: 5 ms

Number of lines: 1

→ ID1: 0 | ID2: 0

Random rate: 0.000000e+00

Normalization factor: 1.000000e+00

ACF: 1.000000e+00

Scatter rate: 0.000000e+00

TOF measurement: 0.000000e+00 ps

TOF measurement range: -1.000000e+00 ps

POI1 (x ; y ; z ) = 0.000000e+00 ; 0.000000e+00 ; -1.000000e+00 ;

POI2 (x ; y ; z ) = 0.000000e+00 ; 0.000000e+00 ; -1.000000e+00 ;
Any help is much appreciated!

[tmerlin]

Something must be wrong in your id calculation if you get an out of range error.
You should check/print the computed IDs in the loop as they are computed and try to understand why you get incorrect values. Most probably at least one of the variable (ringID, rsectorID, nb_crystals, etc…) is wrong.

[ASK]

I’ve checked the cdh file and it seems like only the first event ID has that error with the ID being out of range where as the the rest of the event ID’s are simply 0. This issue only happened when i erased part of the datafileConversion code as mentioned above. When i re-enter those parts that i erased i am just left with a Cdh file with all event ID having ID1 = 0 and ID2 = 0. Is there anyway i can get an example datafileConversion so that i know what i would need to write in the .cc code?

[tmerlin]

Hi,
Originally the IDs in this example have values as the root file

The only other example is the Gate converter code (castor-GATERootToCastor.cc in toolkits), which uses the ConvertIDcylindrical() function implemented in src/management/gDataConversionUtilities.cc. It can be useful to check how it works if you have a Gate version of your system, but otherwise the code may not help you as it is quite generic in order to adapt to the different possible configurations of geometries generated with Gate.
Actually erasing the first part of the loop on nEvents as you did was the right step. You just need to replace them with the ID calculation corresponding to the channel IDs of your original list-mode datafile and to the way you translated your geometry into CASToR (using either the geom file for which case you would have to follow CASToR ID indexing as described in part 5.1.1 of the documentation - or - a LUT build according the 5.1.2 section, which has the advantage that you would not need to follow CASToR indexing as you specify the ID of each crystal).

[ASK]

Sorry, but could you explain in more detail how I would go about writing this conversion? Should it be written In the datafileConversion toolkits code? Im confused since in sec. 5.1.1 the conversion seems to already be written in the gDataConversionUtilities code? Also, i am unsure about needs to be erased in the datafileConversion. As you mentioned the loop of nEvents which i believe is this correct? // Loop on the LORs contained in the event
for(int l=0 ; l<nb_lines_in_event ; l++)
{
// Recover the CASToR ID position of the two crystals according to
// the LUT corresponding to this scanner
castor_id1[l] = ringID1nb_crystals_per_ring
+ rSectorID1nb_crystals_trs
+ crystalID1%nb_crystals_trs ;

castor_id2[l] = ringID2nb_crystals_per_ring
+ rSectorID2nb_crystals_trs
+ crystalID2%nb_crystals_trs ;
}

Is there anything else that would need to be erased in the code or just this one line and everything else stays? Also could you explain the RsectorID1,RsectorID2,ModuleID1,ModuleID2,CrystalID1,CrystalID2 ? They seem to be set to 0 and i am assuming is what is causing my cdh and cdf to be blank (i.e when doing castor-datafileExplore i am getting for every event index ID1:0 and ID2:0). Could you confirm if that is correct? If so how do i fix this issue?

[ASK]

Sorry for all the questions but is there any way to implement a custom LUT into a caster geometry? Based on the general document, I was unsure, but I think it is possible. I have a LUT of our scanner as a csv file but I am not sure how to convert it into the geometry that CASTOR takes, any help is appreciated.

Thank you

[Maxime_Toussaint]

I am not sure what you mean by “implement a custom LUT into a caster geometry”. I will assume that you mean converting your csv file into a CASToR .lut file. In that case, yes it is possible! As shown in p.21 and 22 of CASToR general documentation, a .lut file is simply a binary file of the detector’s position and orientation. You will also need to create a .hscan file (a plain text file) which is used as the header of the .lut file. Look at the examples in config/scanner to learn its format.

[ASK]

Hey maxime,

Yes, I saw that in the general documentation but I’m not sure how to do that given what is said in the general documentation. Is there any way you can guide me on the process? Is there a command where I can specify the CSV file that I would then convert to the binary lut or do I need to write some code?

[Maxime_Toussaint]

You will need to write some code.

In python, I would do something like:
import numpy as np
data = np.genfromtxt(INSERT_HERE_THE_GOOD_OPTIONS)
~Modify data such that it is a 2D array ((N, 6) or (6,N), can’t remember) of detectors position followed by detectors orientation with type float
data.tofile(WHERE_TO_SAVE_THE_DATA)

[tmerlin]

Hi,
The gDataConversionUtilities files only include code to convert data from GATE simulations.
Apart from that, there is no code in CASToR to automatically convert a datafile from a specific system to CASToR format, as each system could have different datafile format (the only exceptions are data converters compiled binaries for several systems, but there is no source code available).
The 5.1.1 section only indicates how castor computes its IDs if you use a geom file. If you want to write your LUT as your last messages suggest, you can forget about this section.

The code in datafileConversion uses ROOT library as example.
In summary, the line Coincidences->GetEntry(i) will update all the values of the root variables linked to Coincidences (RsectorID, ModuleID, CrystalID, etc…), so for each event i the resulting castor_ids will be different according to these variables.
In your case, none of these root variables have any purposes. The GetEntry() function does nothing (actually it is skipped as CASTOR_ROOT is disabled), so you should get rid of all these variables and write your own code instead, specific to your list-mode file, to compute the correct values for castor_ids.

Note that this code is just an example, if it is too confusing or if you prefer to work in another language than c++, you can write your own code from scratch. The only restrictions are the event format described in section 6.1.2.

[ASK]

I have made a LUT and hscan file from my pre-computed LUT (xyz and orientation unit vectors). What steps do I now need to take in order for the datafileConversion? Do I just need to modify the lines

castor_id1[l] = ringID1nb_crystals_per_ring
+ rSectorID1nb_crystals_trs
+ crystalID1%nb_crystals_trs ;

  castor_id2[l] = ringID2*nb_crystals_per_ring
                + rSectorID2*nb_crystals_trs
                + crystalID2%nb_crystals_trs ;

[ASK]

My LUT is formatted in
Scanner element center location (x,y,z): -1.100252e+02 ; -1.263593e+02 ; 8.400000e+01 .Orientation (x,y,z): 6.494480e-01 ; 7.604060e-01 ; 0.000000e+00
as seen from using castor-scannerExplorer with a total of 6144 elements, is that correct and is there any advice on how I would code in the ID of the crystals in the datafileConversion?

Furthermore, for the list mode data, my interpretation of the format from the general documentation is that we need it in the format of 9385.775517368, 5583, 2554 where the first entry is Time of detection in picoseconds (for the event in the trigger region of higher ID number). Next is the channel ID, for the event in the triggered region of higher ID numbers, and then the channel ID for the event in the triggered region of lower ID numbers. And then I have 14,000,000 more coincidences. Would this be the correct format for the list mode data?

Im sorry for the long message but any help is really appreciated.

[tmerlin]

Hi,
Now that you have a LUT, you indeed need to write the function to compute the castor IDs from the data in your original list-mode file.
The number of elements must match the number of discrete detectors in your system. The figures in your LUT seems coherent, but only you can tell if they are consistent with your system.
Regarding the list-mode format, the first entry is time in ms for the detected event. It is mostly used for dynamic reconstruction purposes and is not related to TOF.

[ASK]

Is there any general method or documentation that can guide me in writing the functions to compute the castor IDs from the data in my list-mode file?