In the scenario of low-dose CT image reconstruction, supervised deep learning is widely adopted, which generally demands a dataset containing pairs of normal-dose and correspoding low-dose images and is therefore challenging in clinical situations. This paper proposed an unsupervised deep learning method to tackle this problem.

ref: A Dataset-free Deep Learning Method for Low-Dose CT Image Reconstruction

• Formulation
• Approach
  • Related Works
  • In This Paper
    • Brief Introduction to Variational Approximation
• Method
  • Training
  • Testing
• Summary
  • Why Self-supervised?

Formulation

The image reconstruction problem can be formulated as an inverse problem

where

• $\boldsymbol{A}$ denotes the projection matrix of CT imaging
• $\boldsymbol{y}$ denotes the available measurement
• $\boldsymbol{x}$ denotes the image to be reconstructed
• $n$ denotes the measurement noise, often modeled by i.i.d. random variables

Approach

Related Works

2 approaches to extend unsupervised denoising network to solve inverse problems

1. Treat inverse problems as a denoising process, which post-process reconstructed images
2. Use DIP, which typically implies that early stopping can be an effective technique for regularizing a denoising network

In This Paper

Generally speaking, the proposed method is based on Bayesian inference, where the prior distribution of an image is re-parametrized by a DNN with random weights.

In Bayesian inference, there are 2 representative estimators

1. maximum a posterior (MAP) estimator \(\boldsymbol{x}_{\text{MAP}}=\arg\max\limits_{\boldsymbol{x}} p(\boldsymbol{x}|\boldsymbol{y}).\)
2. MMSE estimator, a.k.a. conditional mean estimatorr \(\boldsymbol{x}_{\text{CM}}=\mathbb{E}_{(\boldsymbol{x}|\boldsymbol{y})}(\boldsymbol{x}|\boldsymbol{y})=\int \boldsymbol{x}\cdot p(\boldsymbol{x}|\boldsymbol{y})\mathrm{d}\boldsymbol{x}.\)

where

denotes the posterior distribution of $\boldsymbol{x}$ given the measurement $y$. It’s overt that both estimators require the posterior distribution, which models the data well.

A common practice in Bayesian inference is to repress it by Bayesian rule

where the likelihood term writes as

in the presence of i.i.d. Gaussian noise $n\sim \mathcal{N}(0, \sigma^2I)$. Then, the study of estimators turn to defining a prior distribution $p(x)$ that accurately describes statistical characters for images to be reconstructed. In this work, the variable $x$ is re-expressed by a DNN with random weights

where $\boldsymbol{x_0}$ is some initial seed and $\boldsymbol{\theta}$ are random variables (netowrk weights). After re-parametrization, the variables for inference now turn to $\boldsymbol{\theta}$; and again, the key to Bayesian inference now is to define an appropriate posterior distribution

for $\boldsymbol{\theta}$. Given that it’s not computationally tractable, variational approximation methods are adopted to approximate the posterior by a set of approximation distributions

parametrized by $\boldsymbol{\mu}$. The optimal approximation with distribution parameters $\boldsymbol{\mu}^{\ast}$ is then estimated by minimizing Kullback–Leibler (KL) divergence.

Brief Introduction to Variational Approximation

The goal of variational estimation is to approximate a posterior, say

with a second distribution, called approximating distribution, $q(\boldsymbol{\beta})$. To make this approximation as accurate as possible, we search over the space of approximating distributions to find the particular distribution with the minimum KL divergence with the actual posterior. Formally, we have to minimize

If no assumptions are made about the factorized distribution, then it is minimized when

which is not helpful, since the posterior itself is generally intractable. So in paper An Introduction to Bayesian Inference via Variational Approximations, the author utilized factorized approximation, because the independence assumption results in the approximating distribution being divided into a set of factors (or blocks of parameters). $\boldsymbol{\beta}$ is first partitioned into a set of $K$ blocks,

Then the approximating distributions turn into

After that, the variational algorithm will identify the specific parametric families that constitute each component of the factorized distribution, which I won’t elaborate on here.

Method

A self-supervised method for LDCT reconstruction from noisy measurement, which is built on the DNN-based re-parametrization for Bayesian inference.

Recall that CT reconstruction can be formulated as the following inverse peoblem:

💡 Given an observed image $\boldsymbol{y}\in \mathbb{R}^m$, which is corrupted according to forward model and noise $n$, find the unkown image $\boldsymbol{x}\in \mathbb{R}^n$ that satisfies the observation

Considering a DNN with random weights for re-parametrization:

Then, the inference of $\boldsymbol{x}$ from noisy measurement $y$ turns to inferring the network weights $\boldsymbol{\theta}$ from $\boldsymbol{y}$. And now we need to derive the posterior, which can be approximated by a set of distributions

where $\odot$ refers to element-wise multiplication, i.e.,

where $\mu_i$ denotes the distribution paramter of $\theta_i$, and $b_i\sim \mathbf{B}(p_i)$ follows a Bernoulli distribution with probability $p_i$, writes as

In other words, the DNN with random weights used in this paper is the widely used network with dropout.

Training

Since we use

to approximate

the optimal estimation is approximated by minimizing the KL divergence between them, i.e.,

For the first term, suppose that $p(\theta)$ is a uniform distribution over a sufficiently large area $\Omega$. For $\theta_i\in {0, 1\times \mu_i}$, we have

and $p(\theta_i)=\frac{1}{s_i}$, where $s_i$ denotes the length of the domain of definition about $\theta_i$. Then the first term tranforms to

It’s apparent from the last derivation that the expression of the first term does not include parameter $\mu$, i.e., it’s a constant.

For the second term, suppose that the measurement noise $n$ is Gaussian so that

then we have

Considering that the first term bears no relevance to parameter $\mu$, we conclude that

It can be seen from the formula above that, the KL divergence only constrains the estimation in the range space¹ of the projection matrix $\boldsymbol{A}$. Recall Eq. (11) and definitions of distribution $q$ and $B(p)$, we deduce that

Integrating total variation (TV) regularization², the loss function for training is defined as

where $\alpha$ is a hyperparameter. Note that when training, dropout is adopted.

Testing

Once the DNN is trained, we have an approximation to the target posterior, denoted by

In this approach, the image $\boldsymbol{x}$ is estimated using the conditional mean estimator, i.e., given measurement $y$, its conditional mean estimator for $\boldsymbol{x}$ reads

By re-parametrization $x=f(\boldsymbol{x_0};\theta)$, we have

By approximating, we have a proximate conditional mean estimator of $\boldsymbol{x}$

This integration is calculated using MC integration ³, i.e., after the network is trained, we take $K$ random samples of the networks with dropout:

Then, the estimate is defined by averaging these $K$ samples

Summary

Why Self-supervised?

To understand why it is formulated in a self-supervised manner, we may start with the plain supervised learning paradigm.

In order to learn for reconstructing images $\boldsymbol{\hat{x}}$ from noisy observations $\boldsymbol{y}$, we generally need a ground truth $\boldsymbol{x}$ to supervise the learning process, i.e., we need to measure the error between the model’s prediction and ground truth to guide the back propagation. For a large amount of data and their corresponding labels ${y_i, x_i}$

However, in real-world scenarios, such perfectly-labeled dataset is hard to acquire. Typically, for medical image reconstruction, if we want supervised learning methodology, labels themselves need to be reconstructed using one algorithm or another, which not only is time-consuming, but also limits the performance and interpretability of the model. Therefore, how to rely merely on the noisy measurements $\boldsymbol{y}$ to train seems both challenging and meaningful. We can turn the formula above into a self-supervised manner by multiplying the projection matrix $A$.

which resembles the formula derived in part Training.