Welcome back to deep learning. So today we want to look into the applications of known operator learning and a particular one that I want to show today is CT reconstruction.

So here you see the formal solution to the CT reconstruction problem and this is the so called filtered back projection or Radon inverse.

This is exactly the equation that I referred to earlier that has already been solved in 1917.

But as you may know CT scanners have only been realized in 1971 so actually Radon who found this very nice solution has never seen it put to practice.

So how did he solve the CT reconstruction problem? Well the CT reconstruction is a projection process and it's essentially a linear system of equations that can be solved.

And the solution is essentially described by a convolution and a sum. So it's a convolution along the detector direction S and then a back projection over the rotation angle theta.

During the whole process we suppress negative values so we kind of also get a non-linearity into the system.

This all can also be expressed in matrix notation so we know that the projection operations can simply be described as a matrix A that describes how the rays intersect with the volume.

And with this matrix you can simply take the volume X multiplied with A and this gives you the projections that you observe in the scanner.

Now getting the reconstruction is you take the projections P and you essentially need some kind of inverse or pseudo inverse of A in order to compute this.

And we can see that there is a solution that is very similar to what we've seen in above continuous equation.

So we have essentially a pseudo inverse here and that is A transpose times A A transpose inverted times P.

And now you could argue that the inverse that you see here in A is actually the filter.

So for this particular problem we know that the inverse of A A transpose will form a convolution.

This is nice because we know how to implement convolutions into deep networks.

So this is what we did. We can map everything into a neural network.

We start on the left hand side. We put in the sinogram and all of the projections.

We have a convolutional layer that is computing the filtered projections.

Then we have a back projection that is a fully connected layer and it's essentially this large matrix A.

And then we have the non-negativity constraint.

So essentially we can define a neural network that does exactly filter back projection.

Now this is actually not so super interesting because there's nothing to learn.

We know all of those weights. And by the way, the matrix A is really huge.

So for 3D problems it can approach up to 65,000 terabytes of memory in floating point position.

So you don't want to instantiate this matrix. And the reason why you don't want to do that is it's very sparse.

So only a very small fraction of the elements in A are actually connections.

So this is very nice for CT reconstruction because then you typically never instantiate A,

but you compute A and A transpose simply using ray tracers.

So this is typically done on the graphics board.

Now why are we talking about all of this?

Well, we've seen there are cases where CT reconstruction is insufficient and we could essentially do trainable CT reconstruction.

But already if you look at the CT book, you already run into the first problems

because if you implement by the book and you just want to reconstruct a cylinder

that is merely showing the value of 1 within this round area,

then you would like to have an image like this one where everything is 1 within the cylinder and outside of the cylinder it's 0.

So we're showing this line plot here along the blue line through the original slice image.

Now if you just implement filtered back projection as you find it in the textbook, you get a reconstruction like this one.

And a typical mistake is that you choose the length of the Fourier transform too short.

And the other one is that you don't consider the discretization appropriately.

Now you can work with this still and fix the problem in the discretization.

So what you can do now is we can essentially train the correct filter.

So what you would do in a classical CT class is you would run through all the math from the continuous integral

to the discrete version in order to figure out the correct filter coefficients.

Instead here we showed that by knowing that it takes the form of a convolution,

we can express our inverse simply as P times the Fourier transform, which is also just a matrix multiplication, right?

Then K is a diagonal matrix that holds the spectral weights.

And then an inverse Fourier transform that is denoted as F Hermitian here.

And then you back project.

So we can simply write this up as a set of matrices.

And by the way, this would then also define as a network architecture.