Dear Castor developers
I ran into this release a week ago and it is really a nice contribution you have made. I have been working with Monte Carlo simulations for many years developing the program simind and often I get this question of supplying reconstruction software since simind only produce projections. So being able to direct to your webpage would be natural.
However, I dont understand how to read data into castor-recon. I also use Interfile and float values so when reading in the manual at you are using interfile for both input data and output gave me confidence for a smooth connection. But I dont see how to actual do this. I checked the spect benchmark file but the cdh format is different. So my question is how I can read in for example a 64x64 float matrices for 64 angles stored in a data file and with an associated interfile header file.
Dear Michael,
Thanks for your interest. Interfile format is only used as input/output format for image (not projection data). As you have seen, input data for SPECT (and also for PET) is composed of 2 files: a header file (.Cdh) and a binary file (.Cdf) storing the data from a machine or a Monte Carlo simulation. Here is some information about this Cdf file. This file is not a pure histogram (a successive number of counts by bin). The Data should (must) be arrange as described in the “CASToR general documentation” page 30.
For instance, in the case of a 64x64 float matrix for 64 angles (and I consider you don’t have both normalization and scatter informations), your Cdf file should contain 262 144 events (64x64x64) with the following format:
EVENT 0:
time: 0 (uint32_t)
data: 32 (float/double) it’s your data for a specific bin
projection ID: 0 (uint32_t) it’s the angle ID between 0 and 63 in this example
binID: 0 (uint32_t) it’s the bin ID between 0 and (64x64 – 1) in this example
We don’t provide a SPECT/PET data converter, but the easiest way to create a correct Cdf file (I consider you are using C/C++) is to create a struct like this:
struct SPECT_CASTOR_EVENT_FORMAT {
uint32_t time;
float data;
uint32_t projectionID;
uint32_t binID;
};
and store each bin from your data in the CASToR event format.
Kind Regards,
Didier Benoit
Dear Didier
Thank you for your explanation. I think it will not be a major problem to implement. For SPECT reconstructions, does the Time value be clearly defined?
Another related question regards the files for the system. I wonder if it would be possible to allow castor-recon program to first look in the local folder for a geometry file before looking in the standard SYSTEM location? If I implement direct writing to the caster format in the simind program, I can foresee that the SYSTEM folder can be quickly filled with files. Alternatively, is it possible to add a path to the system name in the cdh file?
With a time value set to 0, the SPECT reconstruction should work.
In your second question, what do you exactly mean by “SYSTEM location” ? I’m not sure to understand. Is it your scanner system location defined by “CASTOR_CONFIG” during the compilation ? If it’s the case you can’t add a path in the Cdf file for the tag " Scanner name". If you want to create a custom scanner file you have to save it in the “CASTOR_CONFIG” directory.
Kind Regards,
Didier
In our software, we didn’t plan to read the .geom files in 2 or more folders. I think you have 2 solutions :
1- You can copy all your .geom files for each different kind of scanners in the CASTOR_CONFIG folder
2- Or creating a “castor_test” folder, where you write all your .geom files (and the .geom files provided by CASToR) and set CASTOR_CONFIG to the location of “castor_test” during the installation.
Kind regards,
Didier
Hello,
The castor configuration folder also contains other configuration files, so it is not only related to the scanner geometries.
However, a third possibility is to copy the whole 'config' folder (it is light) into your working directory, put your geom files inside the 'scanner' sub-folder and then you can use the '-conf ./your_config_folder' option to overload the path of the configuration directory at run-time.
The '-conf' option is described when calling 'castor-recon -help-misc'.
Regards
Simon
Dear Simon and Didier
Thank you for your comments and information. I have now written a conversion program between simind interfiles and castor and I can reconstruct images 
. I have a final question for the moment and that concerns the units of the mu map. I assume it is mu (cm-1) for the current photon energy. When I looked at the SPECT benchmark file the average value in a ROI is however, around 0.128 cm-1 instead of 0.154 which is this the linear attenuation coefficient for water. Is this an ‘effective attenuation coefficient map’ with lower values with the purpose of compensate for presence of scatter? I just want to be clear about what values the mu-map should reflect.
Dear Michael,
The attenuation coefficient is in cm-1 as explicitly specified in the manual. The attenuation map was extracted directly from the Siemens computation which uses (as you assume) a wide beam approach. The correct value should be the one you mention for the methodology used by CASToR. I apologize for this misunderstanding and a correct mu-map will be uploaded soon.
Thanks for pointing out this small error.
Best regards
Thomas Carlier
Thank you for the clarification Thomas
I still get some problems when trying to reconstruct SIMIND simulated data. Here is an example of a cbf simulation. I have simulated 128 projections of 128x128 and filled the header files as far as I understand. I have a jpeg of the result that seems to include a lot of ring artefacts. Any clues.
Please let me know if someone wants to help me directly without going through the mailing list.
cbf.cdh (1.42 KB)
Thank you for the clarification Thomas
I still get some problems when trying to reconstruct SIMIND simulated data. Here is an example of a cbf simulation. I have simulated 128 projections of 128x128 and filled the header files as far as I understand. I have a jpeg of the result that seems to include a lot of ring artefacts. Any clues.
Please let me know if someone wants to help me directly without going through the mailing list.
OK
It seems to solve the problem but is this a bug in castor not to distribute the angles properly? I use either 64 or 128 projection around 360 degree rotation and always 8 angles/subset so this should be OK.
Another problem I see is if I simulate the cbf with 64 projection around 360. I get very strange images as compared to 128 angles. I have added two jpeg images showing the effect and the related cdh files. Both is reconstructed with –it 10:1
cbf64.cdh (697 Bytes)
cbf128.cdh (1.06 KB)
Dear Michael,
you should pay attention when using OSEM type algorithms. As mentioned in the documentation (page 7) :
The same loop is used whatever the data type (list-mode or histogram) and modality. Thus, when using subsets in the optimizer, it is the user’s responsibility to order the histogram events in an appropriate manner.To avoid potential artefacts in the reconstructed image, we strongly discourage to write the histogram events with the inner loop over the radial or axial coordinate. We recommend to randomly mix the order to the histogram events before writing them in the datafile or, at least, to write the events with the inner loop over the azimuthal coordinate.
Hope this helps,
Best,
Claude Comtat
OK thanks for the information. I did not understand that fully. I will try and rewrite my converter.
Another isssue I found was that reconstructing 128x128 with 128 projections when fine with mlem (-it 10:1) but a 128x128 with 64 projection of exactly the same source and 360 rotation generated very strangle images (note – no osem).
Dear Michael,
radial : I
axial : J
azimuthal : K
You should not store your projection data as it is usually done :
for (k=0 ;k<K ;k++) for (j=0 ;j<J ;j++) for (i=0 ;i<I ;i++) fwrite(Mat(i,j,k));
and then perform an OSEM reconstruction. In CASToR, there is no assumption at all about the organization of the data. The code will read sequentially the data : if there are N subsets, first subset will contain Mat(i=0,j=0,k=0), Mat(i=N,j=0,k=0), Mat (i=2*N,j=0,k=0), … and the resulting subsets will not be well balanced.
You can rather write your matrix as :
for (j=0 ;j<J ;j++) for (i=0 ;i<I ;i++) for (k=0 ;k<K ;k++) fwrite(Mat(i,j,k)) ;
Then, first subset will contain Mat(i=0,j=0,k=0), Mat(i=0,j=0,k=N), Mat(i=0,j=0,k=2*N), … and the resulting subsets will be better balanced.
I did not get your figure 9.
Best,
Claude
Dear Michael,
radial : I
axial : J
azimuthal : K
You should not store your projection data as it is usually done :
for (k=0 ;k<K ;k++) for (j=0 ;j<J ;j++) for (i=0 ;i<I ;i++) fwrite(Mat(i,j,k));
and then perform an OSEM reconstruction. In CASToR, there is no assumption at all about the organization of the data. The code will read sequentially the data : if there are N subsets, first subset will contain Mat(i=0,j=0,k=0), Mat(i=N,j=0,k=0), Mat (i=2*N,j=0,k=0), … and the resulting subsets will not be well balanced.
You can rather write your matrix as :
for (j=0 ;j<J ;j++) for (i=0 ;i<I ;i++) for (k=0 ;k<K ;k++) fwrite(Mat(i,j,k)) ;
Then, first subset will contain Mat(i=0,j=0,k=0), Mat(i=0,j=0,k=N), Mat(i=0,j=0,k=2*N), … and the resulting subsets will be better balanced.
I did not get your figure 9.
Best,
Claude
Thanks
I will try this
I referred to Figure 9 in the general Castor manual – sorry.
I had a similar problem the first time I used CASToR to reconstruct PET data with OSEM.
I’m not familiar with SPECT reconstruction, but I was referring to what is written in the manual, pages 6-7 :
If I follow your convention :
I have a float matrix Mat(I,J,K) where I,J is the projections bins for a specific SPECT camera position (projection ID ) K. The index I goes from left to right and J goes axially from the patient head and down.
As the order of the histogram “events” is important for OSEM, I would recommend the following:
ievent=0;
for (J=1; J<= MtsizJ; J++) {
for (I=1; I<= MtsizI; I++) {
for (K=1; K<= NbViews; K++) {
fwrite(Itime, Mat(I,J,K), K-1,(I-1)+(J-1)*MtxsizI);
ievent++;
}
}
}
instead of the usual:
ievent=0;
for (K=1; K<= NbViews; K++) {
for (J=1; J<= MtsizJ; J++) {
for (I=1; I<= MtsizI; I++) {
fwrite(Itime, Mat(I,J,K), K-1,(I-1)+(J-1)*MtxsizI);
ievent++;
}
}
}
where NbViews is a multiple of the number of subsets.
Hi Claude and others
I finally got it working for osem. Thank you for all help. I have noticed that when using a relatively low angles/subset there appear artefacts in the images that is constant in position though out the number of slices. The example is for 8 angles/subsets. I will investigate effect this a little bit more. Is there any one that has seen this before?
Hälsningar / Best regards
Michael
Dear Michael,
Would it be possible to send us your projections and your log file for this reconstruction?
We will investigate the problem with your data because we did not experience such artifacts.
Thanks
Thomas
