A method to estimate scattered radiation contained in x-ray projections for computed tomography, CT, reconstruction, comprising: constructing an object model (215, 320, 420, 450) based on a plurality of projection images generated by CT scanning of an object using an x-ray radiation source and a detector panel, wherein the object model (215, 320, 420, 450) contains a plurality of voxels, and the constructing of the object model (215, 320, 420, 450) comprises assigning each of the plurality of voxels to a specific material type with density based on each of the plurality of voxels' Hounsfield Units, HU, value; constructing a virtual radiation source (410, 440) based on the x-ray radiation source; constructing a virtual detector panel (430, 460, 520) based on the detector panel; performing a simulated CT scanning of the object model (215, 320, 420, 450) by simulating macroscopic behavior of particles being emitted from the virtual radiation source (410, 440), passing through the object model (215, 320, 420, 450), and being detected by the virtual detector panel (430, 460, 520); and generating a simulated scatter image (245, 630) based on a first subset of particles scattered during the simulated CT scanning of the object model (215, 320, 420, 450); wherein the generating of the simulated scatter image (245, 630) comprises: constructing a Boltzmann Transport Equation, BTE, for the first subset of particles scattered during the simulated CT scanning of the object model (215, 320, 420, 450); using a deterministic method solver to solve the BTE and calculate the corresponding particle fluence distribution values in the plurality of voxels; and generating the simulated scatter image (245, 630) based on the corresponding particle fluence distribution values in the plurality of voxels; characterized in that: the BTE is numerically solved by discretizing phase-space comprising angle, spatial location and/or energy; and the deterministic method solver converges to a solution by refining angular, spatial or energy mesh resolution. constructing an object model (215, 320, 420, 450) based on a plurality of projection images generated by CT scanning of an object using an x-ray radiation source and a detector panel, wherein the object model (215, 320, 420, 450) contains a plurality of voxels, and the constructing of the object model (215, 320, 420, 450) comprises assigning each of the plurality of voxels to a specific material type with density based on each of the plurality of voxels' Hounsfield Units, HU, value; constructing a virtual radiation source (410, 440) based on the x-ray radiation source; constructing a virtual detector panel (430, 460, 520) based on the detector panel; performing a simulated CT scanning of the object model (215, 320, 420, 450) by simulating macroscopic behavior of particles being emitted from the virtual radiation source (410, 440), passing through the object model (215, 320, 420, 450), and being detected by the virtual detector panel (430, 460, 520); and generating a simulated scatter image (245, 630) based on a first subset of particles scattered during the simulated CT scanning of the object model (215, 320, 420, 450); wherein the generating of the simulated scatter image (245, 630) comprises: constructing a Boltzmann Transport Equation, BTE, for the first subset of particles scattered during the simulated CT scanning of the object model (215, 320, 420, 450); using a deterministic method solver to solve the BTE and calculate the corresponding particle fluence distribution values in the plurality of voxels; and generating the simulated scatter image (245, 630) based on the corresponding particle fluence distribution values in the plurality of voxels; characterized in that: the BTE is numerically solved by discretizing phase-space comprising angle, spatial location and/or energy; and the deterministic method solver converges to a solution by refining angular, spatial or energy mesh resolution. constructing a Boltzmann Transport Equation, BTE, for the first subset of particles scattered during the simulated CT scanning of the object model (215, 320, 420, 450); using a deterministic method solver to solve the BTE and calculate the corresponding particle fluence distribution values in the plurality of voxels; and generating the simulated scatter image (245, 630) based on the corresponding particle fluence distribution values in the plurality of voxels; characterized in that: the BTE is numerically solved by discretizing phase-space comprising angle, spatial location and/or energy; and the deterministic method solver converges to a solution by refining angular, spatial or energy mesh resolution. the BTE is numerically solved by discretizing phase-space comprising angle, spatial location and/or energy; and the deterministic method solver converges to a solution by refining angular, spatial or energy mesh resolution.