Iterative reconstruction methods outperform analytical methods when dealing with incomplete or noisy projection data. In iterative reconstruction, forward projection and backprojection are two major operations, the implementation of which depends on the system matrix. Currently, the most used forward projection is Siddon’s method due to its simplicity. By contrast, the Köhler projection can obtain higher precision. However, performing the iterative reconstruction using Köhler’s projection is still a problem to be solved. In this paper, we propose a method for exact computation of the system matrix for Köhler’s projection. First, the intersection points of a ray with a series of grid planes are calculated and merged into an ordered sequence. Then the contribution of each voxel within a cubic gird between two adjacent intersection points are calculated using linear interpolation and Simpson integration. Depending on the adjacent relation between the current and previous cubic grid, we can determine the final weight factors, or accumulate previous voxel contribution into current cubic grid. Experimental results demonstrate that the projection calculated by our method keeps the same precision as that by conventional Köhler’s projection, and the iterative reconstruction quality using Köhler’s projection is significantly improved compared with that using Siddon’s projection.