Transcription
Hello everybody, and welcome to the third and final part in this series looking at imagery construction in CT. We started with simple back projection, moved on to filtered back projection, and now we're going to end off by looking at a final process known as iterative reconstruction.
It's the most commonly used process in modern CT scanners, and this diagram is showing you what we're going to go through today: the steps we're going to take to generate an image, an output image that's of better quality than filtered back projection and requires a lower dose, more importantly.
Now, remember, all of these image reconstruction methods are mathematical processes. We're using projection data to calculate a final image, and in iterative reconstruction, we're going to take input data, enter that data into what's known as an iterative loop, and continually update that input data until we reach our output value.
Now, this process is computationally demanding and it takes longer time than filtered back projection, so we need a reason to do it over filtered back projection. Now, if you're your mind back to filtered back projection, I said there were a couple of shortfalls with FBP. The first is what I focused on the most, which is noise. Filtered back projection doesn't deal well with inherent noise within our projection data, whether that's quantum noise from the photons or the detector noise. It also doesn't account for the system geometry that's actually involved in acquiring the projection data: the fan beam geometry, the size of the focal spot, all of which we're going to look at today. And I've said it doesn't deal well with artifacts. We're going to focus these on these specifically in a future talk where we're looking at artifacts in isolation, but here we can see a ring artifact and beam hardening artifact, which are much more prone to occur in filtered back projection over iterative reconstruction.
Now, I want to take you through some of these processes in filtered back projection, and then when we go through iterative reconstruction, you'll see how that accounts for these shortfalls in filtered back projection. So, if we look at how we go about acquiring the data, the projection data for a CT image using filtered back projection, we need to make some false assumptions in order to create that filtered back projection image. When we use filtered back projection, we are creating what's called a closed-form solution. We've got projection data, we take that projection data, and we have a final output where we've got linear attenuation coefficient values for each and every one of the pixels within our image, and we apply Hounsfield units to those linear attenuation coefficient values.
Now, in order to do we need to make these false assumptions, and we're not accounting for these during filtered back projection. The first assumption that we make is that X-rays originate from a point source, and we know when we look at the focal spot, we rotate this heating machine around here, we can see that the focal spot has some width. All the X-rays aren't coming from some infinitely small point. Now, why is that important? Well, in filtered back projection, if we were to assume X-rays were coming from a point source, they would cast a shadow on our detector that was some form of magnification from our patient here, and see how crisp that shadow is because it's coming from a point source.
Now, we know in reality the focal spot has some width, and we know that X-rays would head from the focal spot towards our patient and cast a shadow. The problem comes is that X-rays are released in an isotropic manner. They're released in all 360 degrees from the focal spot. So, X-rays coming from this point on the focal spot are going to cast a shadow on our detector based on the patient here. The same is going to happen from X-rays coming from the other side of the focal spot, and all the variations in between. What we've got here is what's called a penumbra blurring on the edges of the shadow that's cast on the detector here. This is what's known as geometric blurring.
Now, because we get this geometric blurring based on the fact that X-rays aren't coming from a point source, we go and take that data and in filtered back projection, we propagate that data back over the slice. Remember, we're averaging data back over the slice. As we smear that data out, we're not accounting for the penumbra or the geometric blurring that's caused in the image, and you'll see in iterative reconstruction, we are going to account for this geometry here.
The second assumption that we make when we create a filtered back projection image is that the line integral, that's the linear attenuation coefficients that are summed together in line with one detector, all of those linear attenuation coefficients, we assume that X-rays are passing that through that in a pencil beam geometry, parallel beams to one another. We know what we've got here is fan beam geometry. A single detector is receiving X-rays from a fan beam here. Notice how that beam gets wider and wider as it reaches the detector, so we're not sampling equal voxel widths here. As we head through the patient, that width is getting larger and larger. Again, in filtered back projection, we assume that each voxel is an infinitely small point along this line integral here, and we add up those linear attenuation coefficient values when we go about generating our image. We don't account for the changing geometry.
The third assumption that we make is that the differences in signal that we measure on our detectors in filtered back projection, we say is only due to differences in tissue density. We know that's not true. We know that there is noise generated in the image. For one, that's going to account for random variations in signal intensity, whether it be detector noise or quantum noise coming from the photon interactions with the patient here. If we were to look at noise that we looked at before in filtered back projection, there's one way of reducing that noise to do that in filtered back projection: we can actually increase the dose. We can increase the number of photons that are heading towards the patient, and we do that by increasing the filament current in the cathode. Increasing the filament current in the cathode is going to cause more thermionic emission, more electrons available at the surface of the filament that we can accelerate towards the anode and ultimately produce a larger quantity of X-rays. We're going to increase the quantity of the X-ray beam. The average energy of that beam is the same, but we've got many more X-ray photons heading towards our patient. That's going to increase the dose. With that increase in X-ray photons, we're going to get an increase in signal and an improvement in our signal-to-noise ratio. Watch how that happens here. I'm going to increase dose linearly here. See how the noise is reduced out of our sinogram here, and see how that changes the signal-to-noise ratio in the image that we ultimately produce here.
Now, obviously, we're creating a better image, we're reducing the amount of noise relative to signal, but that comes at the expense of dose. And we know with the ALARA principle, if we can avoid increasing patient dose, especially in say pediatric populations or if we're imaging smaller body parts, if we can prevent this increase in dose, then we need to do that, and that's what iterative reconstruction allows us to do.
Now, you may be wondering, you see the noise in our original data here, can't we just have some algorithm to remove the noise from the projection data before we perform the filtered back projection? And that's a good idea. The issue comes in when we use low doses like we're using here. Noise distribution doesn't follow normal distribution curves like a Gaussian curve. It follows what's known as Poisson noise distribution, where we get skewing of noise or the intensity of noise at lower doses. We get a really poor signal-to-noise ratio at lower doses here. As we increase X-ray photon doses, we get a more normal distribution of noise that we could then remove from this projection data much more easily, and it's this that's the crux of iterative reconstruction, where we go about figuring out attenuation data whilst accounting for this noise distribution, Poisson noise distribution.
Now, the way I like to think of Poisson noise versus Gaussian noise, if you're in a house like mine that has a tin roof, actually over my bedroom, I've got a tree that's over the tin roof, and when it rains very lightly, rain hits the tree, and some of the rain falls through the tree. And if I'm listening to the raindrops hitting the roof at very low doses, that very light rain, that noise is quite random, it's not uniform. We haven't got that low hum. If we were having loads of rain coming down, it's very random noise. And if I had to place a bucket out there and I was catching the rain, it'd be very hard for me to predict how much rain I would catch depending on where I place that bucket. When it's pouring with rain, that noise, that on top of the roof is very much more symmetrical, it's very much more predictable. We've got this continuous noise going on the roof, and if I was to place a bucket, no matter where I place it on the roof, going to catch about the same amount of rain. It's a good way of thinking about it. It's not a perfectly accurate analogy, but it shows you that at these lower doses, noise distribution is much more difficult to calculate, and we actually don't have a closed-form solution like we have in filtered back projection or like we have with Gaussian distribution of noise at higher doses, and we need to use an iterative method where we get closer and closer to the right answer, and that's what I'm going to take you through today.
Now, there are two more false assumptions that we make in filtered back projection, all of which are going to be addressed by iterative reconstruction. The fourth is that it assumes the beam is monoenergetic. We know that the X-ray beam that's heading towards the patient is polyenergetic. It has many different photon energies. There's a range or a spectrum of energies, and if we place anything between our anode here and the detectors here, we're going to get preferential attenuation of the lower energy X-rays. I shouldn't say preferential, no one's choosing to attenuate these lower energy X-rays, it's just more likely that these X-rays will be attenuated due to the photoelectric effect. So, as X-ray photons are passing through a greater and greater thickness of patient, we're getting more and more lower energy X-ray photons being attenuated relative to higher energy X-ray photons. So, as the beam passes through the patient, the average photon energy of the beam gets higher and higher. The beam becomes hardened. It's called beam hardening here, and that's responsible for one of the artifacts that we're going to look at in a later talk. In iterative reconstruction, we can make some compensation for the beam hardening that occurs as the X-rays pass through a patient.
And the final incorrect or false assumption that we make during filtered back projection imaging is we assume that there's an equal exposure throughout the image. If I was to plot the area or the region of the patient that was being imaged by this X-ray beam here, so you've got a Y-axis here and an X-axis here, this is spatial locations, and then I were to rotate the CT machine around that patient, we're going to get differential or differing exposure based on the location within the CT machine. These central regions are going to have a higher relative exposure than the peripheral regions here, and in iterative reconstruction, we can account for this differing exposure here, this exposure variation. This is what's known as a sensitivity image here.
So, now that we've gone through the problems that have occurred in filtered back projection, now let's look at the process of iterative reconstruction and how that goes about addressing those problems. How do we ultimately get a better image in iterative reconstruction CT? Well, the first thing we need to do is input some data into our system. The first piece of data that we're going to acquire is what's known as our measured data. Now, remember, the CT machine doesn't know what the patient anatomy is within it. It only detects X-ray photons that are being transmitted through the patient. So, we rotate around the patient, we gather data, we digitize that data, and place it in a sinogram. Remember, each line in the sinogram here is a different projection, representing a different angle around the patient here, digitized here, digitized detection data, and we call that measured data. That's M, measured data for I, meaning every detector I, is the I detector. Each one of those detectors is going to measure differing photon intensities based on the path that the X-rays have traveled, and these X-ray intensities are independent of one another. They don't relate to one another. It's only related to the line integral, to the attenuation that's occurred along that specific line here. So, this data never changes in iterative reconstruction. We keep our Mi, or our measured data, the same. We've measured that once with the CT machine. This data is going to contain true attenuation data, but it's also going to contain noise, and at low doses, that noise is distributed in Poisson distribution.
The second input that we require in this iterative reconstruction system is what's known as an estimated image. We can place any image in this initial estimate here. Here, I've created a random image that I know is going to vaguely correspond to the image that we're trying to reconstruct, but this could just be a blank image. It could be an empty array, an empty data set that matches the dimensions of the image that we're trying to create. If we look at this more closely, what we could do is theorize: if we were to place this image within our CT machine, what projection data would we get? This is called forward projecting. And now this comes an important part with iterative reconstruction: is we can model the data based on the specific CT machine that we're using. We know the width of the focal spot, we know the geometry of the fan beam, we know the differing X-ray attenuation as it passes through different lengths of our estimated image here. This is calculated mathematically. What we're doing is modeling what the forward projection would be. We're also not getting noise here. We don't have actual detectors, we don't have actual X-ray photons passing through our image here. What we've done is mathematically modeled it. So, the X-ray machine or the CT machine is not actually rotating around our image, but it's a good way to think about how we're creating that data. What we're actually doing is applying what's known as a system matrix, and we're going to mention that by the letter A here. The system matrix accounts for all of those geometry inconsistencies that we saw in filtered back projection, and it's specific to a specific manufacturer's what they want to include in this matrix, as well as a specific CT machine that we're using. As the CT machine changes its parameters, we need to change this matrix. We can think of that matrix going around and creating or generating a sinogram here, a sinogram that's free from noise. We don't have any noise in this image, and it accounts for some of those geometric problems that we were looking at in filtered back projection.
Now, we're not actually rotating a matrix around this. What we're actually doing is taking an array, our matrix, our system matrix in array, and we are multiplying it to this vector image that we've created here, and we're getting what's known as our predicted data set. Here's our predicted data set, and we write that as A of X to the I detector, this predicted data set over all of the detectors here, and that's basically a sum of our matrix applied to our vector. It's our matrix-vector multiplication here over all the detectors, over all the voxels that contribute to signal in one detector through our estimate image here, and we're going to add up all of those predicted values through all the voxels in the image. Maybe a little bit out of the scope of the talk, but you'll see this come up later. So, we've created a second sinogram that is a predicted sinogram based on the estimated image that we've created.
Now, what we need to do is take this input data and somehow apply that data into what's known as an iterative loop, and we want to improve on that initial image estimate and try and closely match it to what type of patient anatomy would ultimately produce this projection data here that we've measured in our CT machine. And to do that, what we're going to do is get a ratio between our measured data and the predicted data that we our system matrix to. We want to see where do these two images differ from one another. The ratio here is going to provide us that difference. The closer these values are together, the closer this number is going to become to one, and the further they're apart, the further away from one it's going to get. We're creating a ratio of differences here.
Now, to do that, we actually apply these data sets over one another. We take these data sets, they match in size, and we divide the measured data by our predicted data here, and we've got a new sinogram. This sinogram represents the ratio of differences between these two data sets. Now, we can take that and enter this data into our iterative loop. What we want to do is take this new data and essentially create a spatial image again in X and Y coordinates as to where the differences are between these two images. We could back project effectively and create a back projected image of these differences here.
Now, instead of back projecting, what we can do is what's known as a transpose of this image. We're essentially mathematically back projecting this sinogram but accounting for those system geometries. We are transposing through that initial matrix that we have applied to our initial estimated image. What that's going to give us is a spatial domain image that represents a essentially back projected image of the ratio between our measured and our predicted sinograms.
Now, I said to you before that depending on spatial location, we're going to get differing exposures because of the rotation of the CT machine around the patient, and we can create what's known as a sensitivity image, and we want to take out those differing exposures from this image that we've created. We only want to see true differences, not differences based on exposure. This is what's known as the sensitivity image. Now, to create the sensitivity image, we can basically fill an array, let's say the number one in all of the array, and then create a transpose or back projection of that array, and that will represent differences purely based on the system matrix that we've chosen here and the varying exposure based on location.
What we've created here by taking this image and subtracting this image from it is what's known as a gradient image. This gradient image is essentially showing us where are the largest differences between our measured data set and our predicted data set based on spatial locations. If you look closely here, you can see that underlying image here. This is what's known as a gradient image, and the formula for this gradient image is this formula here: our ratio between measured and predicted data transposed or back projected minus our sensitivity image. We can then take this gradient image and add it to our initial image estimate. Once we add it to our initial image estimate, that's going to update the input that we are supplying the iterative loop, and I call it X + 1. XK was our first input, now XK + 1 is the second iterative input that we're going to place into the input part of this iterative reconstruction loop. Remember where that's come from: this gradient image that we've generated by creating a ratio between the sinograms, transposing that ratio, taking away the sensitivity image, ultimately creating our gradient image, and we're adding that to our initial image estimate.
What that does then is allows us to update our input here. We can again pass this new updated iteration through the system matrix, AX, now what it's A K+1, our new iteration over all the detectors, and we can repeat this process here, updating the iterative loop, and each time we're seeing what are the differences between the two data sets, we're going to account for the sensitivity image, and then update the new image by those differences.
Now, I know for some of you using the visual explanations here, you can kind of see that this new image is getting closer to what we predict our initial image to be, but it's very difficult one to see how we're accounting for noise in here. I haven't seen any step that's taking out the noise in the image. Surely, if we just start fitting the data closer and closer to one another, we're going to include the noise data here. And secondly, how do we know that this mathematics is sound? What is the basis for say, taking away the sensitivity image here? What's the basis for adding the gradient image onto our initial estimate? And very rarely do I actually go through the mathematics, but in this case, I think it's helpful to look at how we start thinking about modeling the data, how we account for, and how we ultimately get to this formula here where we're adding a gradient image onto our initial estimate, and hopefully this will make sense. Don't worry if the math gets a little bit confusing, especially if you haven't done some higher grade math, but it's the gist that I want to get here. We're going to see the end point of this equation. It's going to lead to the same iterative loop that we looked at first.
What we've got here is what's known as a likelihood formula. It's the Poisson likelihood formula, and it describes the likelihood of our estimate, our input image here, the likelihood that that actually matches the projection data that we're sampling here. What this equation mentions, this part of the equation here, models how Poisson noise distribution occurs in an image when we are comparing the two data sets. This is just a standard formula for the Poisson noise distribution in an image, and this takes the product of those likelihoods and gives us a likelihood figure, a likelihood ratio, as to whether our data actually matches the data that we've measured in the CT machine. The whole goal of iterative reconstruction is to find an image prediction that has the highest likelihood or the highest probability of matching our projected data that we've measured with the CT machine. I hope that makes sense. This just adds up the product of those likelihoods over all of our detector elements.
Now, this is computationally quite difficult. We've got products here, we've got exponents, and what we normally do is we take this Poisson likelihood distribution and we take a log of both sides. We convert this formula into what's known as the log likelihood formula. Now, you may be wondering what this term is here, the projection data factorial. This is just a normalization factor that ensures that these probabilities add up to one. When we take the log of both sides, we get rid of the exponent here, and we get our projection data log, our predicted data, minus our predicted data over all of the detectors, and we get minus some constant here, which is a function of this Mi factorial. We get a log of the likelihood on this side. That's why this is called the log likelihood ratio, and it's no longer a product of all those likelihoods, but it's a sum of all those likelihoods over all the detectors.
Now, because this is a constant and it's not reliant on X itself, we can exclude that from the equation, and this is the log likelihood equation that you might see. Again, don't get caught up. The reason we got log of our projected data here, this is a function of Poisson noise distribution. This is how we're accounting for noise within our image. Again, our goal here is to maximize the likelihood of these two data sets matching. So, what we can do is we can actually plot this function on a graph here. We can say our predicted data, that's the data that we've estimated or put into the system, what's the likelihood that that matches our projection data that we've measured in the CT machine, and we can plot it here.
Now, if we look at this graph, it's very easy to see where that maximum is, and this is the formula for finding that maximum. We want to say what's the X value that will provide the maximum sum of the log likelihoods here, the maximum of this equation, and we know that that point is going to be this point on our graph where the gradient of that graph equals zero. And that's actually the clue here. How we go about finding where the maximum is is by using the gradient of this log likelihood graph. Remember, I said this is not a closed-form solution. We don't have a finite answer. These are likelihoods of probabilities. We need to iterate as we get closer and closer to this maximum log likelihood estimate here, an MLE, you may have heard it before.
So, what we can do is find the gradient of this log likelihood graph here, and we know that if we want to find the gradient of the graph, we take a partial derivative of the function itself. I know this is getting complex, but what we're doing is we're finding the gradient of this graph here. For those of you who actually want to work out this partial derivative, here's the formula for A of X, our system matrix, multiplied by our X vector here, where we've got the system matrix for every detector through all the voxels on that line integral, multiplied by the X vector over all of the voxels in the image. This is the partial derivative here, where we get our system matrix multiplied by the ratio between our measured data and predicted data. You might be thinking, wait, I can see that here. That's what we did initially, minus our projection matrix multiplied by one. You'll see why that becomes important again. We can plot this on our graph here, where the gradient is positive, this plot is going to be positive. Where the gradient is negative, this plot is going to be negative. What we've essentially done here is created our gradient image. We can write this in different terms. We can say the transpose of the ratio between our measured data and predicted data, transpose here, minus a sensitivity image, minus the transpose of one. We've got an invisible one here. This graph here is representing the exact same thing as a gradient image that we created earlier.
So, if we were to enter some image in and say it fell here on the predicted data, remember this is for all the different voxels within the image, we can say that that X value at this point here in the predicted data has a positive gradient. What we want to do is then add that gradient onto our initial estimate, and that's going to move us along this graph here. It's going to move us closer to this maximum likelihood estimate here. If the gradient was of the predicted data was here, we'd see the gradient is negative. If we were to add a negative onto our initial estimate, add a negative, it's going to bring us closer to this MLE. You can see how adding this gradient, adding this gradient image to our initial estimate is going to bring that predicted estimate, that predicted data closer to the maximum likelihood estimate on our log likelihood graph here. For some of you, I may have lost you there. I get that this is complex. It might be a video that you want to go over and revisit some of the steps.
In iterative reconstruction, we account for noise by using the Poisson noise distribution, the likelihood ratio between our predicted and our measured data that we saw in that initial equation here. We take the log of both sides to make it computationally more effective, and then we can find the gradient of this log graph, the gradient image here, the same as what we were doing in this process, and add that to the initial estimate. That's going to bring us closer and closer to the true image that would give us this projection data once the noise has been accounted for. We're not going to overfit the data. These formulas prevent us from overfitting the data and introducing noise into our estimated image because we're applying a system matrix that accounts for that geometry, we're also not going to get that geometric blur. We're also not going to get partial volume averaging because of the fan beam that heads out towards that detector. All of that is accounted for in the system matrix here. Each update on the iterative loop updates that initial estimate and gets us closer and closer to the true anatomy that's causing attenuation in our measured data.
Now, either we go through that cycle a set number of times, say we go through 100 iterations and then spit out an output image, or once the difference, the likelihood difference here, is small enough between our measured data and our estimated data, we can reach a threshold where we're going to output the final image that we're actually going to see on our computer.
So, I hope that made sense to you. We're tying off now, image reconstruction in CT. Remember, they're all mathematical processes. Next, we're going to move on to image dose, the most commonly asked question when it comes to exams. The examiners love to ask about dose because it's very much clinically important. Don't want to expose patients to undue radiation.
Now, if you're studying for a specific exam and you want to prep for those questions, I've linked a question bank in the description box below. Go and check out that question bank. We also cover X-ray, MRI, ultrasound physics. Otherwise, I'll see you all in that next talk where we're going to look at dose in CT imaging. Until then, goodbye everybody.