Serious bug on parallel beam projections

[AnderBiguri]

Actual Behavior

Stripe artefacts appear in the recon

I suspect it is some numerical issue. I recommend using a large DSD/DSO for now, or slightly offseting your angles.

Code to reproduce the problem (If applicable)

% VARIABLE                                   DESCRIPTION                    UNITS
%-------------------------------------------------------------------------------------
% Distances
% This works:
geo.DSD = 2000;                              % Distance Source Detector      (mm)
geo.DSO = 1800;                               % Distance Source Origin        (mm)
% This generate striped images:
geo.DSD = 700;                              % Distance Source Detector      (mm)
geo.DSO = 500; 

geo.nDetector=[640;200];					% number of pixels     400,60         (px)
geo.dDetector=[0.25; 0.25]; 					% size of each pixel            (mm)
geo.sDetector=geo.nDetector.*geo.dDetector; % total size of the detector    (mm)

geo.nVoxel=[640;640;200];                   % number of voxels    400,400,60          (vx)
geo.sVoxel=[160;160;50];                      % total size of the image       (mm)
geo.dVoxel=geo.sVoxel./geo.nVoxel;          % size of each voxel            (mm)

geo.mode='parallel';                        % 'cone' or 'parallel'. Geometry type.


head=headPhantom(geo.nVoxel); 
angles = [0:2:358]*pi/180;

% The following also removes the stripes
% angles=angles+0.01;

proj=Ax(head,geo,angles);

img=OS_SART(proj,geo,angles,1);
plotImg(img,'dim',3)

[tsadakane]

The first and last lines are added.

geo=defaultGeometry('mode', 'parallel');  

% VARIABLE                                   DESCRIPTION                    UNITS
%-------------------------------------------------------------------------------------
% Distances
% This works:
geo.DSD = 2000;                              % Distance Source Detector      (mm)
geo.DSO = 1800;                               % Distance Source Origin        (mm)
% This generate striped images:
geo.DSD = 700;                              % Distance Source Detector      (mm)
geo.DSO = 500; 

geo.nDetector=[640;200];					% number of pixels     400,60         (px)
geo.dDetector=[0.25; 0.25]; 					% size of each pixel            (mm)
geo.sDetector=geo.nDetector.*geo.dDetector; % total size of the detector    (mm)

geo.nVoxel=[640;640;200];                   % number of voxels    400,400,60          (vx)
geo.sVoxel=[160;160;50];                      % total size of the image       (mm)
geo.dVoxel=geo.sVoxel./geo.nVoxel;          % size of each voxel            (mm)

geo.mode='parallel';                        % 'cone' or 'parallel'. Geometry type.

head=headPhantom(geo.nVoxel); 
angles = [0:2:358]*pi/180;

% The following also removes the stripes
% angles=angles+0.01;

proj=Ax(head,geo,angles);

img=OS_SART(proj,geo,angles,1);
plotImg(img,'dim',3)
imwrite(img(:,:,70), 'result.png');

I ran this program several times and found that the result images (img(;,;,70)) are different randomly. The stripe did not always appear.

1. appered..

2. not appeared

3. screen shot of diff (svn diff) of 1 and 2

Envs

• Windows 10
• Matalb R2021a
• Visual Studio 2019
• CUDA 11.2
• GeForce GTX 1060 6GB

[tsadakane]

FYI. So far, On python, the following code does not reproduce this phenomena on my PC:

import numpy as np
import matplotlib.pyplot as plt
import tigre
import tigre.algorithms as algs
from tigre.utilities import sample_loader
from tigre.utilities.Measure_Quality import Measure_Quality
import tigre.utilities.gpu as gpu

geo = tigre.geometry(mode="parallel", nVoxel=np.array([200, 640, 640]), default=True)

# VARIABLE                                   DESCRIPTION                    UNITS
#-------------------------------------------------------------------------------------
# Distances
# This works:
geo.DSD = 2000                              # Distance Source Detector      (mm)
geo.DSO = 1800                              # Distance Source Origin        (mm)
# This generate striped images:
geo.DSD = 700                               # Distance Source Detector      (mm)
geo.DSO = 500 

geo.nDetector=np.array([200, 640])          # number of pixels     400,60         (px)
geo.dDetector=np.array([0.25, 0.25])        # size of each pixel            (mm)
geo.sDetector=geo.nDetector * geo.dDetector # total size of the detector    (mm)

geo.nVoxel=np.array([200, 640, 640])        # number of voxels    400,400,60          (vx)
geo.sVoxel=np.array([50, 160, 160])         # total size of the image       (mm)
geo.dVoxel=geo.sVoxel / geo.nVoxel          # size of each voxel            (mm)

geo.mode="parallel"                         # 'cone' or 'parallel'. Geometry type.

head=sample_loader.load_head_phantom(geo.nVoxel)
# angles = [0:2:358]*pi/180
angles = np.arange(start=0, stop = 360, step = 2, dtype=np.float32)*np.pi/180

# The following also removes the stripes
# angles=angles+0.01

proj=tigre.Ax(head,geo,angles)

img=algs.ossart(proj, geo, angles, 1)
plt.imsave("result.png", img[70])
tigre.plotImg(img, dim="Z")

[AnderBiguri]

Thanks @tsadakane for testing!

I suspect this may have to do with a missing memset, due to its erratic behaviour. Its likely that here a cudamemset is missing and that adding it would fix the problem. Will try when I have access to my computer:

https://github.com/CERN/TIGRE/blob/master/Common/CUDA/Siddon_projection_parallel.cu#L281

EDIT: Or maybe somewhere else.... This does not seem like a place where a memset would be needed.

[AnderBiguri]

retry, I may have fixed this by fixing #453