Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

BUFGIX: correct TOF kernel truncation (taking tof bin width into account) #70

Merged
merged 14 commits into from
Jun 14, 2024
Merged

BUFGIX: correct TOF kernel truncation (taking tof bin width into account) #70

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
14 commits
Select commit Hold shift + click to select a range
c4df3dc
add tof_norm wip script
gschramm Jun 13, 2024
cb8548a
add sigma eff to C TOF projectors to do better TOF truncation
gschramm Jun 13, 2024
8c213ac
add effective sigmas to CUDA projectors
gschramm Jun 13, 2024
2fb9ec8
add new test to test whether sum over tof is close to nontof
gschramm Jun 13, 2024
1de14ad
change version to 1.8.1
gschramm Jun 13, 2024
83195d4
add global correction for TOF kernel truncation to C projectors
gschramm Jun 13, 2024
9157233
add global TOF truncation factor to CUDA projectors
gschramm Jun 13, 2024
05b5276
use tigher tolerance in tof sum test
gschramm Jun 13, 2024
436c225
print cuda kernel file in __init__
gschramm Jun 14, 2024
fa775e3
add better test to test sino / LM forward projects comparing against …
gschramm Jun 14, 2024
daf7702
BUGFIX: add missing devices in pet_scanner phis + test
gschramm Jun 14, 2024
de1635b
print cuda flags on import
gschramm Jun 14, 2024
7880a91
make sure cuda_kernel_file is defined
gschramm Jun 14, 2024
f553a69
add extra tests that test exception for non-int tofbin arrays
gschramm Jun 14, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
36 changes: 21 additions & 15 deletions c/src/joseph3d_back_tof_lm.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -66,6 +66,17 @@ void joseph3d_back_tof_lm(const float *xstart,
float sig_tof = (lor_dependent_sigma_tof == 1) ? sigma_tof[i] : sigma_tof[0];
float tc_offset = (lor_dependent_tofcenter_offset == 1) ? tofcenter_offset[i] : tofcenter_offset[0];

// calculate the effective sigma (standard deviation of the TOF Gaussian convolved with the tofbin width)
// we need to to device which TOF bins a certain voxel along an LOR
float sig_eff = sqrtf(sig_tof*sig_tof + tofbin_width*tofbin_width/12.0f);

// factor that corrects the sum of the TOF weights to be 1, assuming that tofbin_width << sig_tof
// for n_sigma = 3.5, this factor is 1.0004
// for n_sigma = 3, this factor is 1.0027
// for n_sigma = 2.5, this factor is 1.0126
// for n_sigma = 2, this factor is 1.0476
float tof_trunc_corr_factor = 1.0f / erff(n_sigmas/sqrtf(2));

float xstart0 = xstart[i*3 + 0];
float xstart1 = xstart[i*3 + 1];
float xstart2 = xstart[i*3 + 2];
Expand Down Expand Up @@ -178,8 +189,8 @@ void joseph3d_back_tof_lm(const float *xstart,
//-- check where we should start and stop according to the TOF kernel
//-- the tof weights outside +- 3 sigma will be close to 0 so we can
//-- ignore them
istart_tof_f = (x_m0 + (it*tofbin_width - n_sigmas*sig_tof)*u0 - img_origin0) / voxsize0;
iend_tof_f = (x_m0 + (it*tofbin_width + n_sigmas*sig_tof)*u0 - img_origin0) / voxsize0;
istart_tof_f = (x_m0 + (it*tofbin_width - n_sigmas*sig_eff)*u0 - img_origin0) / voxsize0;
iend_tof_f = (x_m0 + (it*tofbin_width + n_sigmas*sig_eff)*u0 - img_origin0) / voxsize0;

if (istart_tof_f > iend_tof_f){
tmp = iend_tof_f;
Expand Down Expand Up @@ -228,8 +239,8 @@ void joseph3d_back_tof_lm(const float *xstart,
powf((x_m1 + (it*tofbin_width + tc_offset)*u1 - x_v1), 2) +
powf((x_m2 + (it*tofbin_width + tc_offset)*u2 - x_v2), 2));

//calculate the TOF weight
tw = 0.5f*(erff((dtof + 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)) -
// calculate the TOF weight including the correction factor due to kernel truncation
tw = 0.5f*tof_trunc_corr_factor*(erff((dtof + 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)) -
erff((dtof - 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)));

if ((i1_floor >= 0) && (i1_floor < n1) && (i2_floor >= 0) && (i2_floor < n2))
Expand Down Expand Up @@ -301,8 +312,8 @@ void joseph3d_back_tof_lm(const float *xstart,
//-- check where we should start and stop according to the TOF kernel
//-- the tof weights outside +- 3 sigma will be close to 0 so we can
//-- ignore them
istart_tof_f = (x_m1 + (it*tofbin_width - n_sigmas*sig_tof)*u1 - img_origin1) / voxsize1;
iend_tof_f = (x_m1 + (it*tofbin_width + n_sigmas*sig_tof)*u1 - img_origin1) / voxsize1;
istart_tof_f = (x_m1 + (it*tofbin_width - n_sigmas*sig_eff)*u1 - img_origin1) / voxsize1;
iend_tof_f = (x_m1 + (it*tofbin_width + n_sigmas*sig_eff)*u1 - img_origin1) / voxsize1;

if (istart_tof_f > iend_tof_f){
tmp = iend_tof_f;
Expand Down Expand Up @@ -352,7 +363,7 @@ void joseph3d_back_tof_lm(const float *xstart,
powf((x_m2 + (it*tofbin_width + tc_offset)*u2 - x_v2), 2));

//calculate the TOF weight
tw = 0.5f*(erff((dtof + 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)) -
tw = 0.5f*tof_trunc_corr_factor*(erff((dtof + 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)) -
erff((dtof - 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)));

if ((i0_floor >= 0) && (i0_floor < n0) && (i2_floor >= 0) && (i2_floor < n2))
Expand Down Expand Up @@ -420,13 +431,8 @@ void joseph3d_back_tof_lm(const float *xstart,
//-- check where we should start and stop according to the TOF kernel
//-- the tof weights outside +- 3 sigma will be close to 0 so we can
//-- ignore them
istart_tof_f = (x_m1 + (it*tofbin_width - n_sigmas*sig_tof)*u1 - img_origin1) / voxsize1;

//-- check where we should start and stop according to the TOF kernel
//-- the tof weights outside +- 3 sigma will be close to 0 so we can
//-- ignore them
istart_tof_f = (x_m2 + (it*tofbin_width - n_sigmas*sig_tof)*u2 - img_origin2) / voxsize2;
iend_tof_f = (x_m2 + (it*tofbin_width + n_sigmas*sig_tof)*u2 - img_origin2) / voxsize2;
istart_tof_f = (x_m2 + (it*tofbin_width - n_sigmas*sig_eff)*u2 - img_origin2) / voxsize2;
iend_tof_f = (x_m2 + (it*tofbin_width + n_sigmas*sig_eff)*u2 - img_origin2) / voxsize2;

if (istart_tof_f > iend_tof_f){
tmp = iend_tof_f;
Expand Down Expand Up @@ -476,7 +482,7 @@ void joseph3d_back_tof_lm(const float *xstart,
powf((x_m2 + (it*tofbin_width + tc_offset)*u2 - x_v2), 2));

//calculate the TOF weight
tw = 0.5f*(erff((dtof + 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)) -
tw = 0.5f*tof_trunc_corr_factor*(erff((dtof + 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)) -
erff((dtof - 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)));

if ((i0_floor >= 0) && (i0_floor < n0) && (i1_floor >= 0) && (i1_floor < n1))
Expand Down
44 changes: 29 additions & 15 deletions c/src/joseph3d_back_tof_sino.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -82,6 +82,17 @@ void joseph3d_back_tof_sino(const float *xstart,
float istart_tof_f, iend_tof_f;
int istart_tof, iend_tof;

// calculate the effective sigma (standard deviation of the TOF Gaussian convolved with the tofbin width)
// we need to to device which TOF bins a certain voxel along an LOR
float sig_eff = sqrtf(sig_tof*sig_tof + tofbin_width*tofbin_width/12.0f);

// factor that corrects the sum of the TOF weights to be 1, assuming that tofbin_width << sig_tof
// for n_sigma = 3.5, this factor is 1.0004
// for n_sigma = 3, this factor is 1.0027
// for n_sigma = 2.5, this factor is 1.0126
// for n_sigma = 2, this factor is 1.0476
float tof_trunc_corr_factor = 1.0f / erff(n_sigmas/sqrtf(2));

// test whether the ray between the two detectors is most parallel
// with the 0, 1, or 2 axis
d0 = xend0 - xstart0;
Expand Down Expand Up @@ -208,13 +219,13 @@ void joseph3d_back_tof_sino(const float *xstart,

// get the relevant tof bins (the TOF bins where the TOF weight is not close to 0)
relevant_tof_bins(x_m0, x_m1, x_m2, x_v0, x_v1, x_v2, u0, u1, u2,
tofbin_width, tc_offset, sig_tof, n_sigmas, n_half,
tofbin_width, tc_offset, sig_eff, n_sigmas, n_half,
&it1, &it2);

for(it = it1; it <= it2; it++){
//--- add extra check to be compatible with behavior of LM projector
istart_tof_f = (x_m0 + (it*tofbin_width - n_sigmas*sig_tof)*u0 - img_origin0) / voxsize0;
iend_tof_f = (x_m0 + (it*tofbin_width + n_sigmas*sig_tof)*u0 - img_origin0) / voxsize0;
istart_tof_f = (x_m0 + (it*tofbin_width - n_sigmas*sig_eff)*u0 - img_origin0) / voxsize0;
iend_tof_f = (x_m0 + (it*tofbin_width + n_sigmas*sig_eff)*u0 - img_origin0) / voxsize0;

if (istart_tof_f > iend_tof_f){
tmp = iend_tof_f;
Expand All @@ -233,8 +244,9 @@ void joseph3d_back_tof_sino(const float *xstart,
powf((x_m1 + (it*tofbin_width + tc_offset)*u1 - x_v1), 2) +
powf((x_m2 + (it*tofbin_width + tc_offset)*u2 - x_v2), 2));

//calculate the TOF weight
tw = 0.5f*(erff((dtof + 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)) -
// calculate the TOF weight
// correct for the fact that the sum of the TOF weights is not 1
tw = 0.5f*tof_trunc_corr_factor*(erff((dtof + 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)) -
erff((dtof - 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)));

if ((i1_floor >= 0) && (i1_floor < n1) && (i2_floor >= 0) && (i2_floor < n2))
Expand Down Expand Up @@ -338,13 +350,13 @@ void joseph3d_back_tof_sino(const float *xstart,

// get the relevant tof bins (the TOF bins where the TOF weight is not close to 0)
relevant_tof_bins(x_m0, x_m1, x_m2, x_v0, x_v1, x_v2, u0, u1, u2,
tofbin_width, tc_offset, sig_tof, n_sigmas, n_half,
tofbin_width, tc_offset, sig_eff, n_sigmas, n_half,
&it1, &it2);

for(it = it1; it <= it2; it++){
//--- add extra check to be compatible with behavior of LM projector
istart_tof_f = (x_m1 + (it*tofbin_width - n_sigmas*sig_tof)*u1 - img_origin1) / voxsize1;
iend_tof_f = (x_m1 + (it*tofbin_width + n_sigmas*sig_tof)*u1 - img_origin1) / voxsize1;
istart_tof_f = (x_m1 + (it*tofbin_width - n_sigmas*sig_eff)*u1 - img_origin1) / voxsize1;
iend_tof_f = (x_m1 + (it*tofbin_width + n_sigmas*sig_eff)*u1 - img_origin1) / voxsize1;

if (istart_tof_f > iend_tof_f){
tmp = iend_tof_f;
Expand All @@ -363,8 +375,9 @@ void joseph3d_back_tof_sino(const float *xstart,
powf((x_m1 + (it*tofbin_width + tc_offset)*u1 - x_v1), 2) +
powf((x_m2 + (it*tofbin_width + tc_offset)*u2 - x_v2), 2));

//calculate the TOF weight
tw = 0.5f*(erff((dtof + 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)) -
// calculate the TOF weight
// correct for the fact that the sum of the TOF weights is not 1
tw = 0.5f*tof_trunc_corr_factor*(erff((dtof + 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)) -
erff((dtof - 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)));

if ((i0_floor >= 0) && (i0_floor < n0) && (i2_floor >= 0) && (i2_floor < n2))
Expand Down Expand Up @@ -468,13 +481,13 @@ void joseph3d_back_tof_sino(const float *xstart,

// get the relevant tof bins (the TOF bins where the TOF weight is not close to 0)
relevant_tof_bins(x_m0, x_m1, x_m2, x_v0, x_v1, x_v2, u0, u1, u2,
tofbin_width, tc_offset, sig_tof, n_sigmas, n_half,
tofbin_width, tc_offset, sig_eff, n_sigmas, n_half,
&it1, &it2);

for(it = it1; it <= it2; it++){
//--- add extra check to be compatible with behavior of LM projector
istart_tof_f = (x_m2 + (it*tofbin_width - n_sigmas*sig_tof)*u2 - img_origin2) / voxsize2;
iend_tof_f = (x_m2 + (it*tofbin_width + n_sigmas*sig_tof)*u2 - img_origin2) / voxsize2;
istart_tof_f = (x_m2 + (it*tofbin_width - n_sigmas*sig_eff)*u2 - img_origin2) / voxsize2;
iend_tof_f = (x_m2 + (it*tofbin_width + n_sigmas*sig_eff)*u2 - img_origin2) / voxsize2;

if (istart_tof_f > iend_tof_f){
tmp = iend_tof_f;
Expand All @@ -493,8 +506,9 @@ void joseph3d_back_tof_sino(const float *xstart,
powf((x_m1 + (it*tofbin_width + tc_offset)*u1 - x_v1), 2) +
powf((x_m2 + (it*tofbin_width + tc_offset)*u2 - x_v2), 2));

//calculate the TOF weight
tw = 0.5f*(erff((dtof + 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)) -
// calculate the TOF weight
// correct for the fact that the sum of the TOF weights is not 1
tw = 0.5f*tof_trunc_corr_factor*(erff((dtof + 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)) -
erff((dtof - 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)));

if ((i0_floor >= 0) && (i0_floor < n0) && (i1_floor >= 0) && (i1_floor < n1))
Expand Down
35 changes: 23 additions & 12 deletions c/src/joseph3d_fwd_tof_lm.c
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -68,6 +68,17 @@ void joseph3d_fwd_tof_lm(const float *xstart,
float sig_tof = (lor_dependent_sigma_tof == 1) ? sigma_tof[i] : sigma_tof[0];
float tc_offset = (lor_dependent_tofcenter_offset == 1) ? tofcenter_offset[i] : tofcenter_offset[0];

// calculate the effective sigma (standard deviation of the TOF Gaussian convolved with the tofbin width)
// we need to to device which TOF bins a certain voxel along an LOR
float sig_eff = sqrtf(sig_tof*sig_tof + tofbin_width*tofbin_width/12.0f);

// factor that corrects the sum of the TOF weights to be 1, assuming that tofbin_width << sig_tof
// for n_sigma = 3.5, this factor is 1.0004
// for n_sigma = 3, this factor is 1.0027
// for n_sigma = 2.5, this factor is 1.0126
// for n_sigma = 2, this factor is 1.0476
float tof_trunc_corr_factor = 1.0f / erff(n_sigmas/sqrtf(2));

float xstart0 = xstart[i*3 + 0];
float xstart1 = xstart[i*3 + 1];
float xstart2 = xstart[i*3 + 2];
Expand Down Expand Up @@ -181,8 +192,8 @@ void joseph3d_fwd_tof_lm(const float *xstart,
//-- check where we should start and stop according to the TOF kernel
//-- the tof weights outside +- 3 sigma will be close to 0 so we can
//-- ignore them
istart_tof_f = (x_m0 + (it*tofbin_width - n_sigmas*sig_tof)*u0 - img_origin0) / voxsize0;
iend_tof_f = (x_m0 + (it*tofbin_width + n_sigmas*sig_tof)*u0 - img_origin0) / voxsize0;
istart_tof_f = (x_m0 + (it*tofbin_width - n_sigmas*sig_eff)*u0 - img_origin0) / voxsize0;
iend_tof_f = (x_m0 + (it*tofbin_width + n_sigmas*sig_eff)*u0 - img_origin0) / voxsize0;

if (istart_tof_f > iend_tof_f){
tmp = iend_tof_f;
Expand Down Expand Up @@ -249,8 +260,8 @@ void joseph3d_fwd_tof_lm(const float *xstart,
powf((x_m1 + (it*tofbin_width + tc_offset)*u1 - x_v1), 2) +
powf((x_m2 + (it*tofbin_width + tc_offset)*u2 - x_v2), 2));

//calculate the TOF weight
tw = 0.5f*(erff((dtof + 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)) -
// calculate the TOF weight, including the correction factor for kernel truncation
tw = 0.5f*tof_trunc_corr_factor*(erff((dtof + 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)) -
erff((dtof - 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)));

p[i] += (tw * cf * toAdd);
Expand Down Expand Up @@ -299,8 +310,8 @@ void joseph3d_fwd_tof_lm(const float *xstart,
//-- check where we should start and stop according to the TOF kernel
//-- the tof weights outside +- 3 sigma will be close to 0 so we can
//-- ignore them
istart_tof_f = (x_m1 + (it*tofbin_width - n_sigmas*sig_tof)*u1 - img_origin1) / voxsize1;
iend_tof_f = (x_m1 + (it*tofbin_width + n_sigmas*sig_tof)*u1 - img_origin1) / voxsize1;
istart_tof_f = (x_m1 + (it*tofbin_width - n_sigmas*sig_eff)*u1 - img_origin1) / voxsize1;
iend_tof_f = (x_m1 + (it*tofbin_width + n_sigmas*sig_eff)*u1 - img_origin1) / voxsize1;

if (istart_tof_f > iend_tof_f){
tmp = iend_tof_f;
Expand Down Expand Up @@ -367,8 +378,8 @@ void joseph3d_fwd_tof_lm(const float *xstart,
powf((x_m1 + (it*tofbin_width + tc_offset)*u1 - x_v1), 2) +
powf((x_m2 + (it*tofbin_width + tc_offset)*u2 - x_v2), 2));

//calculate the TOF weight
tw = 0.5f*(erff((dtof + 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)) -
// calculate the TOF weight including the correction factor for kernel truncation
tw = 0.5f*tof_trunc_corr_factor*(erff((dtof + 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)) -
erff((dtof - 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)));


Expand Down Expand Up @@ -418,8 +429,8 @@ void joseph3d_fwd_tof_lm(const float *xstart,
//-- check where we should start and stop according to the TOF kernel
//-- the tof weights outside +- 3 sigma will be close to 0 so we can
//-- ignore them
istart_tof_f = (x_m2 + (it*tofbin_width - n_sigmas*sig_tof)*u2 - img_origin2) / voxsize2;
iend_tof_f = (x_m2 + (it*tofbin_width + n_sigmas*sig_tof)*u2 - img_origin2) / voxsize2;
istart_tof_f = (x_m2 + (it*tofbin_width - n_sigmas*sig_eff)*u2 - img_origin2) / voxsize2;
iend_tof_f = (x_m2 + (it*tofbin_width + n_sigmas*sig_eff)*u2 - img_origin2) / voxsize2;

if (istart_tof_f > iend_tof_f){
tmp = iend_tof_f;
Expand Down Expand Up @@ -486,8 +497,8 @@ void joseph3d_fwd_tof_lm(const float *xstart,
powf((x_m1 + (it*tofbin_width + tc_offset)*u1 - x_v1), 2) +
powf((x_m2 + (it*tofbin_width + tc_offset)*u2 - x_v2), 2));

//calculate the TOF weight
tw = 0.5f*(erff((dtof + 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)) -
// calculate the TOF weight including the correction factor for kernel truncation
tw = 0.5f*tof_trunc_corr_factor*(erff((dtof + 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)) -
erff((dtof - 0.5f*tofbin_width)/(sqrtf(2)*sig_tof)));

p[i] += (tw * cf * toAdd);
Expand Down
Loading
Loading