# CT corrupted by Poisson noise, weighted least squares, how to calculate in detail

I've trained a convolutional neural network to recognize ellipses. I took random generated ellipses as ground truth images and applied radon and inverse radon transform for the noisy images.

Now I want to add Poisson noise to the sinogram before applying the inverse radon transform. For the reconstruction, the paper I'm following suggests to replace the data-fidelity term in

$${min}_x \frac{1}{2} || Hx-y||^2_2$$

with weighted least squares. (Here, x is the input image and y is what we are measuring).

It goes on with: "Let $$\mu$$ represent the distribution of linear attenuation coefficient of an object and $$[A \mu]_m$$ represents their line integral. The mth CT measurement, $$y_m$$ is a Poisson random variable with parameters

$$p_m \sim Poisson( b_m e^{-[A \mu]_m} +r_m)$$

$$y_m = -log(\frac{p_m}{b_m})$$

where $$b_m$$ is the blank scan factor and $$r_m$$ is the readout noise."

The paper then suggests to use these $$p_m$$ to build a diagonal weight matrix $$W$$ to get

$${min}_x \frac{1}{2} || W^{1/2}Hx-W^{1/2}y||^2_2$$.

But what exactly is $$[A \mu]_m$$? Where do I get $$A$$ and $$\mu$$? In a linked paper, it says

$$[A \mu]_m=\sum_{j=1}^J a_{mj} \mu_j, m=1, \dots , M$$, where $$A=\{a_{ij}\}$$ is the system matrix and M and J is the the number of projections and pixels.

But again, what exactly is the system Matrix $$A$$, where can I find it?

Here are the papers I'm referring to: number 1 (page 12 (1451), part B) and number 2 (page 3 (259), part B).

Any help is highly appreciated, thank you very much.

Edit:

Thanks to @A_A, I now have an idea of how $$A$$ looks like. But I still don't really know how to calcalute $$p_m$$. Would it be possible to get $$p_m$$ from $$y_m=−log(\frac{p_m}{b_m})$$, that means $$p_m=b_m e^{−y_m}$$? But how exactly would I determine $$y_m$$; let's say my image $$x$$ is from $$\mathbb{R}^{512×512}$$ and performing a radon tranformation gives me a sinogram $$Hx \in \mathbb{R}^{512×32}$$. Then $$W \in \mathbb{R}^{512×512}$$ and there are $$512$$ values of $$p_m$$. How would I go on from there? And what should I choose for the blank scan factor?

Thank you very much.

• Please see updated response. – A_A Dec 13 '18 at 15:44

The first paper is more clear than the second.

"Let $$\mu$$ represent the distribution of linear attenuation coefficient of an object and $$[A \mu]_m$$ represents their line integral. "

The signal obtained from Computed Tomography (CT) depends on how much was the X-Ray beam attenuated as it travels through the human body (there are other modalities such as the Positron Emission Tomography (PET) where the signal depends on how much emission you measure. Same reconstruction technique, different signal model).

You don't really know what is inside some $$A$$. It represents the body. What you do know is that it comes in the form of a matrix of sites which is called $$A$$. Now that you have the $$A$$, then you consider the total attenuation along the path of a ray through this $$A$$. That is the "line integral" $$[A \mu]_m$$ .

The $$m^{th}$$ CT measurement, $$y_m$$ is a Poisson random variable with parameters $$p_m \sim Poisson( b_m e^{-[A \mu]_m} +r_m)$$ $$y_m = -log(\frac{p_m}{b_m})$$ where $$b_m$$ is the blank scan factor and $$r_m$$ is the readout noise."

This $$[A \mu]_m$$ is inaccessible to you. What is accessible to you is $$y_m$$, which is what you measure, it contains noise and is expressed here as the ratio of a measurement through the body divided by a measurement in air ("blank scan factor").

Notice here, $$y_m$$ is the line integral. One number. One ray through $$A$$. That is, one pixel from one line of your Radon transform. One measurement. NOT the whole line of the Radon transform. To look along the line of the Radon you need to do it in an iteration (i.e, expressions where $$\sum_{m=1}^M y_m \times ...$$).

Hope this helps.

EDIT:

I still don't really know how to calcalute $$p_m$$.

$$p_m$$'s expression is given in the paper as coming from a Poisson distribution. It does not look very much Poisson because it is missing the rate ($$k$$) and so it ends up being closer to an exponentially decaying distribution. The effect this is trying to simulate is the variation of the detector's sensitivity with attenuation (or radiation energy).

Would it be possible to get $$p_m$$ from $$y_m=−log(\frac{p_m}{b_m})$$, that means $$p_m=b_m e^{−y_m}$$?

That means $$p_m \sim Poisson( b_m e^{-[A \mu]_m} + r_m)$$. $$r_m$$ is the readout noise.

But how exactly would I determine $$y_m$$;

[...]

And what should I choose for the blank scan factor?

I don't know what should you choose for $$b_m, r_m$$, it depends on what you are modeling (detector voltage? Hounsfield scale?). Usually, $$A$$ is supposed to contain Hounsfield scale values directly, in which case your blank scan factor should be that for air. But this means that the rest of your image values are adjusted according to what sort of tissue they depict. The point of this is to take into account the Poisson noise because it modulates the readings. So, perhaps you can lump things together into reasonable values obtained by measurements over phantoms.

let's say my image $$x$$ is from $$\mathbb{R}^{512×512}$$ and performing a radon tranformation gives me a sinogram $$Hx \in \mathbb{R}^{512×32}$$.

Just a note here, if you keep your "detector width" at 512 measurements, you will be sampling a disc inscribed within the $$512 \times 512$$ square. When the detector is looking at a $$512 \times 512$$ image at a $$\frac{\pi}{4}$$ bearing, its diagonal is $$\approx 724$$ pixels. What I like to do is to center the square image in a black frame that is $$diagonal \times diagonal$$ in dimensions and use that for reconstruction. It is also easier to keep the detector geometry fixed and rotate the image. For a parallel ray detector, the line integrals are simply the sums of the columns of the image. For a fan beam, each ray deviates but you don't have to re-calculate this for every projection.

Then $$W \in \mathbb{R}^{512×512}$$ and there are $$512$$ values of $$p_m$$. How would I go on from there?

As far as I can tell, this $$W$$ term is $$w_m = \frac{(p_m - r_m)^2}{p_m}$$ (equation 29).

But, the part of the discussion on page 1451 that involves $$W$$ is to explain how to calculate your $$E$$ that drives the optimisation detailed in equation 33. between what you "guess" $$Hx$$ (see "data fidelity term" pages 1440-1441 around equations 1 and 2) and what you measure $$ym$$.

So how image is affected by Poisson noise is one thing and calculating that error term during optimisation is another.

Hope this helps.

• @CarlosDanger Yes, I think that the second paper could do better. I am not sure I follow with "...source and measurement locations and the pixel definitions..." (?). In a CT scanner, the geometry is fixed (?). – A_A Dec 12 '18 at 19:39
• Thank you very much for your detailed answer! So, it's possible to calculate $p_m$ from $y_m=-log(\frac{p_m}{b_m}$, that means $p_m=b_m e^{-y_m}$? But how exactly would I determine $y_m$; let's say my image $x$ is from $\mathbb{R}^{512 \times 512}$ and performing a radon tranformation gives me a sinogram $Hx \in \mathbb{R}^{512 \times 32}$. Then $W \in \mathbb{R}^{512 \times 512}$ and there are 512 values of $p_m$. How would I go on from there? And what should I choose for the blank scan factor? I'm kinda lost and would be grateful for another answer. Thanks. – M. Evans Dec 12 '18 at 19:41
• @M.Evans Can I please ask you to edit your original question with this additional information because it was not there when I provided this answer and I have focused on $A$. – A_A Dec 13 '18 at 11:24
• @CarlosDanger I hope we are not talking cross purposes here :) Yes, I agree, somehow we have to represent the data that comes from the scanner. This is what $A$ does. I am not sure that these two sources use the same notation however. The authors of the cited paper in this question (p 1451), are doing a quadratic approximation of $E$ which already involves $Hx$ that is already the line integral. The line integral would have already taken into account the $w$s that the algebraic reconstruction book is talking about. – A_A Dec 13 '18 at 11:30
• I've edited my question, thank you very much for your help. – M. Evans Dec 13 '18 at 12:54

The two papers use slightly different terminology for $$A$$.

Tomography works by measuring attenuation coefficients through an object (in dB units, so the values add up linearly) by irriadiating it from many different angles and taking measurements on the other side.

In this picture there is only one source, so to get a good image, the source rotates around the circle and takes measurements at every angle. This creates a series of measurements that look something like this: In the first paper, this notation is used: $$\mathbf y = \mathbf H \mathbf x + \mathbf n$$ where $$\mathbf x$$ is a vector of attenuation coefficients, one per pixel $$j$$, in some dB-like units, $$\mathbf y$$ is the measured attenuation along each path $$i$$, and $$\mathbf H$$ is an $$M \times N$$ matrix representing the mapping along each path where you basically add up the attenuation pixels, multiplied by some weight that represents what "fraction" of each pixel is "covered" by that ray. You know how much power the source is emitting, so you obtain $$\mathbf y$$ by taking the received power at each detector location, divide by the emitted power, and then convert to dB.

This equation is called the "forward model", in that it is a model of how the physics of the hidden object $$\mathbf x$$ produce what is observable $$\mathbf y$$.

You basically have to construct $$\mathbf H$$ manually from the geometry by drawing lines between all the source and detector locations and seeing which pixels are crossed. In this formulation, $$\mathbf H$$ and $$\mathbf y$$ are known, and $$\mathbf x$$ is unknown. In this paper, $$\mathbf A$$ represents some linear operator that is used to estimate $$\mathbf x$$, i.e. $$\hat{\mathbf x} = \mathbf A \mathbf y = \ \mathbf A \mathbf H \mathbf x$$ This equation is called the "inverse model" or "inverse problem", in that you are attempting to invert the physics of the system to estimate the vector of hidden parameters $$\mathbf x$$ from the measurements $$\mathbf y$$. Many algorithms exist for linear reconstruction, such as algebraic reconstruction technique (ART), which basically is just using $$\mathbf A$$ to do a least-squares estimate, or filtered back projections (FBP), which uses Fourier methods based on the Radon transform. In the first paper, $$\mathbf A$$ is stated to be an FBP algorithm (below eq 14).

In the second paper, the equation for the forward model is: $$\mathbf p = \textrm A \mu$$ where $$\mathbf p$$ is the vector of measurements, $$\mathrm A$$ is the system matrix, and $$\mu$$ is the vector of attenuation coefficients. This yields to the following correspondence between the two papers:

• (Paper 1) $$\leftrightarrow$$ (Paper 2)
• $$\mathbf y \leftrightarrow \mathbf p$$
• $$\mathbf H \leftrightarrow \mathrm A$$
• $$\mathbf x \leftrightarrow \mu$$