Error in fanbeam strip kernel

[murdock-grewar]

Hi. There appears to be a bug in astra's fan-beam projection code. I have included image examples and source code to produce a particularly potent example.

Examples

First, here is the volume being projected (it is nonnegative everywhere):

Here is a sinogram I produced with since-lost settings. Each row corresponds to a different projection angle. There are discontinuities between rows. Here is another sinogram, where the issue is much more prevalent. There are even negative values, as can be seen from the colour scale on the right-hand side.

Source code

Here is a .py file which reproduces the second sinogram:

#!/usr/bin/env python
import astra
import cupy as cp
import numpy as np
import cupyx.scipy.ndimage
import cupyx.profiler
import scipy.ndimage
import scipy.special
import cupyx.scipy.fft as cufft

from typing import *

import pyqtgraph as pg
from pyqtgraph.Qt import QtWidgets

def plot_attenuation(
      data: Union[np.ndarray,cp.ndarray], 
      init_zslicenum: Optional[int] = None, 
      title: Optional[str] = None,
      show: bool = True,
      ):
    # Interpret image data as row-major instead of col-major
    pg.setConfigOptions(imageAxisOrder='row-major')

    app = pg.mkQApp("Attenuation Plot")

    ## Create window with ImageView widget
    win = QtWidgets.QMainWindow()
    win.resize(1200, 800)
    imv = pg.ImageView(discreteTimeLine=True)
    win.setCentralWidget(imv)
    win.show()
    win.setWindowTitle('pyqtgraph plot' if title is None else title)
    imv.setHistogramLabel("Attenuation")

    if type(data) is np.ndarray:
      dnp = data
    else:
      dnp = cp.asnumpy(data)
    imv.setImage(dnp)
    imv.setLevels(min=np.amin(dnp), max=np.amax(dnp))
    imv.setCurrentIndex(dnp.shape[0]//2 if init_zslicenum is None else init_zslicenum)

    # Start up with an ROI
    imv.ui.roiBtn.setChecked(False)
    imv.roiClicked()

    if show:
      pg.exec()
    return app, win, imv

def square_frame(
        inner_diam_pix: int,
        outer_diam_pix: int,
        imgwid: int,
    ):
    att = np.zeros((imgwid,imgwid), dtype=np.float64)
    empty_wid = (imgwid-outer_diam_pix)//2
    noncenter_wid = (imgwid-inner_diam_pix)//2
    att[
        empty_wid:imgwid-empty_wid,
        empty_wid:imgwid-empty_wid,
    ] = 1
    att[
        noncenter_wid:imgwid-noncenter_wid,
        noncenter_wid:imgwid-noncenter_wid,
    ] = 0
    return att

def phantom(wid: int) -> np.ndarray:
    ph1 = square_frame(
        inner_diam_pix = int(wid * 0.2),
        outer_diam_pix = int(wid * 0.4),
        imgwid = wid,
    ) 
    ph2 = square_frame(
        inner_diam_pix = int(wid * 0.60),
        outer_diam_pix = int(wid * 0.69),
        imgwid = wid,
    ) 
    return ph1 + ph2

def soften_from_rad(relr_fad):
    def soften_fun_r(r_arr: np.ndarray) -> np.ndarray:
        return (1/2)*(1+np.cos(
            np.clip(r_arr - relr_fad, 0, 1 - relr_fad) * (
                (np.pi) / (1 - relr_fad)
            )
        ))
    return soften_fun_r 

def main():

    # Config
    wid:     int = 351
    numangs: int = 1000
    SO           = wid/2*(1.2)#(1.08)
    OD           = wid

    # Calc
    det_wid = 5 + 2*int(np.tan(np.arcsin(wid/2/SO))*(SO+OD))
    fan_angle = 2*np.arcsin(wid/2/SO)
    print(det_wid, fan_angle*180/np.pi)

    # Main
    Theta = np.linspace(start=0, stop=2*np.pi, endpoint=False, num=numangs)
    vol_geom = astra.create_vol_geom(wid, wid)
    lst_proj_geom = [
        astra.create_proj_geom(
            'fanflat', 
            1, 
            det_wid, 
            [theta + np.pi/2], 
            SO, 
            OD,
        ) for theta in Theta
    ]
    lst_projector_ids = [
        astra.create_projector(
            'strip_fanflat', 
            # 'line_fanflat',
            proj_geom,
            vol_geom,
        ) for proj_geom in lst_proj_geom
    ]

    # Funcs
    def A_i(i: int, data: np.ndarray) -> np.ndarray:
        id, sinogram = astra.create_sino(data, lst_projector_ids[i])
        return sinogram
    def A(data: np.ndarray) -> np.ndarray:
        return np.stack(
            [A_i(i, data)/numangs for i in range(numangs)],
            axis = 0,
        )

    x = phantom(wid)
    Ax = A(x)
    plot_attenuation(x[np.newaxis,:,:], title='phantom')
    plot_attenuation(Ax[np.newaxis,:,0,:], title='Ax')

if __name__ == '__main__':
    main()

[murdock-grewar]

I had a quick look at the fan beam projection code.

I believe this issue may have something to do with the switch_t and associates, because the discontinuities are at 3*PIdiv4, 5*Pidiv4, 7*PIdiv4, 9*PIdiv4 and so on.

https://github.com/astra-toolbox/astra-toolbox/blob/7273498d63451ee75c5389b693fe301f9f04d356/include/astra/FanFlatBeamStripKernelProjector2D.inl#L88-L95

[askorikov]

I can reproduce the issue, but not quite sure about the origin yet. line_fanflat projector seems to work fine.