For the $i^\textrm{th}$ ray, $r_i$, the measured attenuation is $\rho_i$:
$\rho_i = \int \mu(t) \; dt$
over the path of $r_i$.
Let $f(x,y)$ by the 2D slice or $f(x,y,z)$ be the 3D body. In either case, $f$ is divided into pixels (2D) or voxels (3D) and these elements are labelled sequentially:
$f_1, f_2, f_3, \ldots, f_N$
for $N$ elements in all.
$f_j$ is the attenuation (measured in absorption per length) in the $j^\textrm{th}$ element.
Here are $N$ elements and a ray, $r_i$:
There are $M$ rays. For ray $r_i$, let $w_{ij}$ be the weight of element $f_j$ with respect to $r_i$. That is, $w_{ij}$ is the relative importance of $f_j$ to the attenuation computation along $r_i$:
$\rho_i = w_{i1} f_1 + w_{i2} f_2 + \ldots + w_{iN} f_N = \sum_j w_{ij} f_j$
If ray $r_i$ does not intersect element $f_j$, the corresponding weight is zero: $w_{ij} = 0$.
If the weight, $w_{ij}$, is the length of the intersection of ray $r_i$ with element $f_j$, then
$w_{ij} \times f_j = $ (length) $\times$ (attenuation per length)
which is the total attenuation along $r_i$ as it passes through $f_j$.
There are $N$ unknowns (the $f_1, \ldots, f_N$) and $M$ equations (one for each of the rays: $r_1, \ldots r_M$):
$W \; f = \rho$
where $W$ is a $M$-row by $N$-column matrix of the $w_{ij}$, $f$ is a $N$-row vector of the $f_j$, and $\rho$ is a $M$-row vector of the $\rho_i$.
$ \left[ \begin{array}{ccccccc} w_{11} & w_{12} & w_{13} & \cdots & w_{1N} \\ w_{21} & w_{22} & & & \vdots \\ w_{31} \\ \vdots \\ w_{M1} & \cdots & & & w_{MN} \end{array} \right] \; \left[ \begin{array}{c} f_1 \\ f_2 \\ \vdots \\ f_N \end{array} \right] \; = \; \left[ \begin{array}{c} \rho_1 \\ \rho_2 \\ \rho_3 \\ \vdots \\ \rho_M \end{array} \right] $
This is a big system of linear equations, and can be solved by many methods.
But it's too big (e.g. $M = 512 \times 512 = 262144$ columns in matrix) to solve by anything other than an iterative method.
Kaczmarz's Method
Consider the problem $W \; f = \rho$ with known $W$ and $\rho$, and unknown $f$.
Let $f^{(k)}$ be the $k^\textrm{th}$ approximation to the solution, $f$.
Kaczmarz's method computes successive approximations $f^{(1)}, f^{(2)}, f^{(3)}, \ldots$.
Each $f^{(k)}$ is a point in an $N$-dimensional space of solutions.
Each row of $W \; f = \rho$ is a plane equation:
$\displaystyle\rho_i = w_1 f_1 + w_2 f_2 + \ldots + w_N f_N = \sum_j w_{ij}\ f_j = w_i \cdot f$
To compute $f^{(k+1)}$ from $f^{(k)}$, the point $f^{(k)}$ is projected onto the "next" plane. Then $f^{(k+1)}$ is projected onto the next plane after that, and so on.
The initial $f$ is typically $0$. The method iterates through one equation row (i.e. plane), then the next equation row (i.e. plane), and so on.
Point Projection
To project a point, $f$, onto a plane $w \cdot f = \rho$:
If we are projecting onto the $i^\textrm{th}$ row:
$f^{(k+1)} = f^{(k)} - { f^{(k)} \cdot w_i - \rho_i \over w_i \cdot w_i } \; w_i $
Above, $\Delta f = { f^{(k)} \cdot {w_i \over |w_i|} - {\rho_i \over |w_i|}}$ is the error distance (along ray $r_i$) between the projection, $f^{(k)} \cdot {w_i \over |w_i|}$, of our current solution and the observed projection, $\rho_i$.
That error, $\Delta f$, is normalized with respect to $|w|$ and is then distributed over the $f$ in proportion to the normalized element weights, $w_i \over |w|$.
Simplifications
Ideally, $w_{ij}$ is the area of intersection of $f_j$ with the beam of non-zero width of ray $r_i$. But length works as well, and the units of $f_j$ in "absorption per length" match when using length.
To simplify (for faster computation), we can approximate $w_{ij}$ as
$w_{ij} \approx \left\{ \begin{array}{ll} 1 & \textrm{if beam crosses "near" centre of } f_j \\ 0 & \textrm{otherwise} \end{array} \right.$
Acceleration
Iterate through rays in random order, since $r_i$ and $r_{i+1}$ are almost the same and the projections onto the corresponding planes do not differ very much.