# How to compute the weight matrix for tomography applications?

I am trying to compute the weight matrix for a set of straight ray paths through a reconstruction region.

Ideally, I would like to be able to do this for both a rectangular grid region, where each grid cell is a rectangular prism, and a geodetic grid region, where grid cell boundaries are specified in latitude, longitude, and altitude.

I have looked into the ASTRA library, but cannot get it to compile on my system.

I have tried implementing the algorithm myself, but it is rather non-robust, buggy, and I feel like there should be a spelt out algorithm somewhere to do this.

## 1 Answer

The canonical way to do this is using Sidons algorithm.

Where you go from there depends on your application. If your ray paths have finite width, you may want to consider that in your reconstruction.

For the geodesics, things get more complicated. I'm fairly confident that an analytical solution exists, but I simply don't know it. I would suggest solving this by brute-force instead. Simply trace your line through your geometry in small steps, check which region each step is in and add the distance between this step and the previous to corresponding matrix index.. This is essentially doing the integration via the trapezoidal method.

I'm doing something similar in the reverse case (curved ray path, regtangular grid) and it works extremely well in practice.

If your data or reconstruction region gets larger than a very modest size (say, 256x256x64), you may want to consider never forming the matrix explicitly, and simply implementing the matrix-vector product directly, for use in iterative reconstruction.