Feature reconstruction from incomplete tomographic data without detour

In this paper, we consider the problem of feature reconstruction from incomplete x-ray CT data. Such problems occurs, e.g., as a result of dose reduction in the context medical imaging. Since image reconstruction from incomplete data is a severely ill-posed problem, the reconstructed images may suffer from characteristic artefacts or missing features, and significantly complicate subsequent image processing tasks (e.g., edge detection or segmentation). In this paper, we introduce a novel framework for the robust reconstruction of convolutional image features directly from CT data, without the need of computing a reconstruction firs. Within our framework we use non-linear (variational) regularization methods that can be adapted to a variety of feature reconstruction tasks and to several limited data situations . In our numerical experiments, we consider several instances of edge reconstructions from angularly undersampled data and show that our approach is able to reliably reconstruct feature maps in this case.

READ FULL TEXT VIEW PDF

page 14

page 15

page 16

page 17

page 18

01/27/2020

Medical image reconstruction with image-adaptive priors learned by use of generative adversarial networks

Medical image reconstruction is typically an ill-posed inverse problem. ...
06/13/2019

Deep Variational Networks with Exponential Weighting for Learning Computed Tomography

Tomographic image reconstruction is relevant for many medical imaging mo...
12/03/2018

JSR-Net: A Deep Network for Joint Spatial-Radon Domain CT Reconstruction from incomplete data

CT image reconstruction from incomplete data, such as sparse views and l...
02/22/2022

Data-Consistent Local Superresolution for Medical Imaging

In this work we propose a new paradigm of iterative model-based reconstr...
12/24/2020

Parallel-beam X-ray CT datasets of apples with internal defects and label balancing for machine learning

We present three parallel-beam tomographic datasets of 94 apples with in...
02/26/2020

3D Phase Retrieval at Nano-Scale via Accelerated Wirtinger Flow

Imaging 3D nano-structures at very high resolution is crucial in a varie...
04/10/2021

Shape reconstructions by using plasmon resonances

We study the shape reconstruction of a dielectric inclusion from the far...

1 Introduction

Computed tomography (CT) has established itself as one of the standard diagnostic tools in medical imaging. However, the relatively high radiation dose that is used to produce high resolution CT images (and that patients are exposed to) has become one of the major clinical concerns [CTDoseReduction, cancer, thousands, middleage]. The reduction of the radiation exposure of a patient while ensuring the diagnostic image quality constitutes one of the main challenges in CT. Beside patient safety, also the reduction of scanning times and costs constitute important aspects of dose reduction, which is often achieved by reducing the x-ray energy level (leading to higher noise levels in the data) or by reducing the number of collected CT data (leading to incomplete data), cf. [CTDoseReduction].

In this work, we particularly consider incomplete data situations, e.g., that arise in a sparse or limited view setup, where CT data is collected only with respect to a small number of x-rays directions or within a small angular range. The intentional reduction of the angular sampling rate leads to an under-determined and severely ill-posed image reconstruction problem, c.f. [natterer]. As a consequence, the reconstructed image quality can be substantially degraded, e.g., by artefacts or missing features [juergentodd], and complicate subsequent image processing tasks (such as edge detection or segmentation) that are often employed within a CAD-pipeline (computer aided diagnosis). Therefore, the development of robust feature detection algorithms for CT that ensure the diagnostic image quality is an important and very challenging task. In this paper, we introduce a framework for feature reconstruction directly from incomplete tomographic data, which is in contrast to the classical 2-step approach where reconstruction and feature detection are performed in two separate steps.

Incomplete tomographic data

In this article, we consider the parallel beam geometry and use the 2D Radon transform as a model for the (full) CT data generation process, where denotes the unit sphere in and is a function representing the sought tomographic image (CT scan). Here, the value represents one x-ray measurement over a line in

that is parametrized by the normal vector

and the signed distance from the origin . In what follows, we consider incomplete data situations where the Radon data are available only for a small number of directions, given by . We denote the angularly sampled tomographic Radon data by . In this context, the (semi-discrete) CT data will be called incomplete if the Radon transform is insufficiently sampled with respect to the directional variable. Prominent instances of incomplete data situations are: sparse angle setup, where the directions in are sparsely distributed over the full angular range ; limited view setup, where covers only small part of the full angular range . Precise mathematical criteria of (in-) sufficient sampling can be derived from the Shannon sampling theory. Those criteria are based on the relation between the number of directions and the bandwidth of , cf. [natterer]

. In this work, we will mainly focus on the sparse angle case, with uniformly distributed directions

on a half-sphere, e.g., directions with uniformly distributed angles .

Feature reconstruction in tomography

In the following, we consider image features that can be extracted from a CT scan by a convolution with a kernel . In this context, the notion of a feature map will refer to the convolution product , and the convolution kernel

will be called feature extraction filter. Examples of feature detection tasks that can be realized by a convolution include edge detection, image restoration, image enhancement, or texture filtering

[jain1989fundamentals]. For example, in case of edge detection, the filter can be chosen as a smooth approximation of differential operators, e.g., of the Laplacian operator [jahne2005digital]. In our practical examples, we will mainly focus on edge detection in tomography. However, the proposed framework also applies to more general feature extraction tasks.

In many standard imaging setups, image reconstruction and feature extraction are realized in two separate steps. However, as pointed out in [louis2], this 2-step approach can lead to unreliable feature maps since feature extraction algorithms have to account for inaccuracies that are present in the reconstruction. This is particularly true for the case of incomplete CT data as those reconstructions may contain artefacts. Hence, combining these two steps into an approach that computes feature maps directly from CT data can lead to a significant performance increase, as was already pointed out in [louis1, louis2]. In this work, we account for this fact and extend the results of [louis1, louis2] to more general setting and, in particular, to limited data situations.

Main contributions and related work

In this work, we propose a framework to directly reconstruct the feature map from the measured tomographic data. Our approach is based on the forward convolution identity for the Radon transform, that is , where on the right hand side the convolution is taken with respect to the second variable of the Radon transform, cf. [natterer]. This identity implies that, given (semi-discrete) CT data, the feature map satisfies the (discretized) equation , where is the modified (preprocessed) CT data. Therefore, the sought feature map can be formally computed by applying a discretized version of the inverse Radon transform to , i.e., as . In case of full data (sufficient sampling), this can be accurately and efficiently computed by using the well known filtered backprojection (FBP) algorithm with the filter . However, if the CT data is incomplete, this approach would lead to unreliable feature maps since in such situations the FBP is known to produce inaccurate reconstruction results, cf. [natterer, juergentodd].

In order to account for data incompleteness, we propose to replace the inverse by a suitable regularization method for that is also able to deal with undersampled data. More concretely, we propose to reconstruct the (discrete) feature map by minimizing the following Tikhonov type functional:

This framework offers a flexible way for incorporating a-priori information about the feature map into the reconstruction and, in this way, to account for the missing data. For example, from the theory of compressed sensing it is well known that sparsity can help to overcome the classical Nyquist-Shannon-Whittaker-Kotelnikov paradigm [cs]. Hence, whenever the sought feature map is known to be sparse (e.g., in case of edge detection), sparse regularization techniques can be easily incorporated into this framework.

Approaches that combine image reconstruction and edge detection have been proposed for the case of full tomographic data, e.g., in [louis1, louis2]. Although the presented work follows the spirit of [louis1, louis2], it comes with several novelties and advantages. On a formal level, our approach is based on the forward convolution identity, in contrast to the dual convolution identity, given by , that is employed in in [louis1, louis2]. The latter requires full (properly sampled) data, since the backprojection operator integrates over the full angular range (requiring proper sampling in the angular variable). In contrast, our framework is applicable to incomplete Radon data situations, since the forward convolutional identity (used in our approach) can be applied to more general situations. Moreover, in order to recover the feature map , we use non-linear regularization methods that can be adapted to a variety of situations and incorporate different kinds of prior information. From this perspective, our approach also offers more flexibility. A similar approach was presented in our recent proceedings article [frikel2021combining], where the main focus was on the stable recovery of the image gradient from CT data and its application to Canny edge detection. Following the ideas of [louis1, louis2] similar feature detection methods were developed also for other types of tomography problems, e.g., in [Hahn_2013, Rigaud_2015, Rigaud2017]. Besides that, we are not aware of any further results concerning convolutional feature reconstruction from incomplete x-ray CT data.

Combinations of reconstruction and segmentation have also been presented in the literature for different types of tomography problems, e.g., in [elangovan2001sinograms, klann2011mumford, storath2015joint, Burger_2016, Romanov2016, LQWJ2018, JoinRecSeg2018]. As a commonality to our approach, many of those methods are based on the minimization of an energy functional of the form , incorporating feature maps as prior information. Important examples include Mumford-Shah like approaches [klann2011mumford, Burger_2016, JoinRecSeg2018, LQWJ2018] or the Potts model [storath2015joint]. Also, geometric approaches for computing segmentation masks directly from tomographic data were employed in [elangovan2001sinograms].

Outline

In Section 2 we provide some basic facts about the Radon transform, sampling and sparse recovery. In Section 3 we introduce the proposed feature reconstruction framework and present several examples of convolutional feature reconstruction filters along with corresponding data filters, mainly focusing on the case of edge detection. Experimental results will be presented in Section 4. We conclude with a summary and outlook given in Section 5.

2 Materials and Methods

In this section, we recall some basic facts about the 2D Radon transform, including important identities and sampling conditions. In particular, we define the sub-sampled Radon transform that will be used throughout this article. Although, our presentation is restricted to the 2D case (because this makes the presentations more concise and clear), the presented concepts can be easily generalized to the -dimensional setup.

2.1 The Radon transform

Let denote the Schwartz space on (space of smooth functions that are rapidly decaying together with all their derivatives) and the Schwartz space over as the space of all smooth functions that are rapidly decaying together with all their derivatives in the second component, cf. [natterer]. We consider the Radon transform as an operator between those Schwartz spaces, , which is defined via

(2.1)

where , and denotes the rotated version of by counterclockwise (in particular, is a unit vector perpendicular to ). The value represents one x-ray measurement along the x-ray path that is given by the line . Since , the following symmetry property holds for the Radon transform, . Hence, it is sufficient to know the values of Radon transform only on a half sphere, e.g., on the upper half sphere. Such data is therefore considered to be complete. The dual transform (backprojection operator) is defined as ,

(2.2)

The Radon transform is a well defined linear and injective operator, and several analytic properties are well-known. One of the most important properties is the so-called Fourier slice theorem that describes the relation between the Radon and the Fourier transforms. In order to state this relation, we first recall that the Fourier transform is defined as

, for . Whenever convenient, we will also use the abbreviated notation . The Fourier transform is a linear isomorphism on the Schwartz space and its inverse is given by . In what follows, we will denote the convolution of two functions by , where . Moreover, for functions , the Fourier transform will refer to the 1D-Fourier transform of with respect to the second variable. Analogeously, will denote the convolution of with respect to the second variable.

Lemma 2.1 (Properties of the Radon transform).

  1. [label=(R0), leftmargin=3em]

  2. Fourier slice theorem: .

  3. Convolution identity: .

  4. Dual convolution identity: .

  5. Intertwining with Derivatives:

  6. Intertwining with Laplacian: .

Proof.

All identities are derived in [natterer, Chapter II]. ∎

The approach that we are going to present in Section 3 is based on the convolution identity and can be formulated for arbitrary spatial dimension . For the sake of clarity we consider two spatial dimensions in this paper. In this case, we will use the parametrization of given by with . Then . For the Radon transform we will (with some abuse of notation) write

2.2 Sampling the Radon transform

Since in practice one has to deal with discrete data, we are forced to work with discretized (sampled) versions of the Radon transform. In this context, questions about proper sampling arise. In order to understand what it means for the CT data to be complete (properly sampled) or incomplete (improperly sampled), we recall some basic facts from the Shannon sampling theory for the Radon transform for the case of parallel scanning geometry (see for example [natterer, Section III]).

In what follows, we assume that is compactly supported on the unit disc and consider sampled CT data with equispaced angles in and equispaced values in for the -variable, i.e.,

(2.3)

For given sampling points (2.3) and a finite dimensional subspace we define the discrete Radon transform as

(2.4)

The basic question of classical sampling theory in the context of CT is to find conditions on the class of images and on the sampling points under which the sampled data uniquely determines the unknown function . Sampling theory for CT has been studied, for example, in [desbat1993efficient, Far04, Far06, natterer95sampling, RatLin81]. While the classical sampling theory (e.g., in the setting of classical signal processing) works with the class of band-limited functions, the sampling conditions in the context of CT are typically derived for the class of essentially band-limited function.

Remark 2.2 (Band-limited and essentially band-limited functions).

A function is called -band-limited if its Fourier transform vanishes for . A function is called essentially -band-limited if is negligible for in the sense that is sufficiently small, see [natterer]. One reason for working with essentially band-limited functions in CT is that functions with compact support cannot be strictly band-limited. However, the quantity can become arbitrarily small for functions with compact support.

The bandwidth is crucial for the correct sampling conditions and the calculation of appropriate filters. If consists of essentially -band-limited functions that vanish outside the unit disc , then the correct sampling conditions are given by [natterer]

(2.5)

Obviously, as the bandwidth increases, the step sizes and have to decrease in order that (2.5) is satisfied. Thus, if the bandwidth is large, a large number measurements (roughly ) have to be collected. As a consequence, for high resolution imaging the sampling conditions require a large number of measurements. Thus, in practical applications, high resolution imaging in CT also leads to large scanning times and to high doses of x-ray exposure. A classical approach for dose reduction consists in reduction of x-ray measurements. For example, this can be achieved by angular undersampling where Radon data is collected only for a relatively small number of directions .

Definition 2.3 (Sub-sampled Radon transform).

Let be defied by (2.5) and let be the set of essentially -band-limited functions that vanishes outside the unit disc (note that in that case the discrete Radon transform defined in (2.4) is correctly sampled). For we call

(2.6)

the sub-sampled discrete Radon transform. We will also use the semi-discrete Radon transform , where we only sample in the angular direction but not in the radial direction.

If we perform actual undersampling, where the number of directions in is much less than , then the linear equation will be is severely under-determined and its solution requires additional prior information (e.g., sparsity of the feature map).

3 Feature reconstruction from incomplete data

In this section, we present our approach for feature map reconstruction from incomplete data. For a given bandwidth , we let denote the set of essentially -band-limited functions that vanishes outside . Furthermore, we assume that the set of directions is chosen according to the sampling conditions (2.5).

Problem 1 (Feature reconstruction task).

Let and let be the noisy subsampled (semi-discrete) CT data with , where is the true but unknown image and the known noise level. Given a feature extraction filter

, our goal is to estimate the feature map

from the (undersampled) data .

Remark 3.1.
  1. [itemsep=1ex]

  2. From a general perspective, the Problem 1 is related to the field of optimal recovery [micchelli1977survey] where the goal is to estimate certain features of an element in a space from noisy indirect observations.

  3. Depending on the particular choice of the filter , the Problem 1 corresponds to several typical tasks in tomography. For example, if is chosen as an approximation of the Delta distribution, the Problem 1 is equivalent to the classical image reconstruction problem. In fact, the filtered backprojection algorithm (FBP) is derived in this way from the dual convolution identity 3 for the full data case, cf. [natterer]. Another instance of Problem 1 is edge reconstruction from tomographic data . For example, this can be achieved by choosing the feature extraction filter as the Laplacian of an approximation to the Delta distribution (e.g., Laplacian of Gaussian (LoG)). Then, the Problem 1 boils down to an approximate recovery of the Laplacian of , which is used in practical edge-detection algorithms (e.g., LoG-filter [jain1989fundamentals, jahne2005digital]).

  4. Traditionally, the solution of Problem 1 is realized via the 2-step approach: First, estimate and, second, apply convolution in order to estimate the feature map . This 2-step approach has several disadvantages: Since image reconstruction in CT is (possibly severely) ill-posed, the fist step might introduce huge errors in the reconstructed image. Those errors will also be propagated through the second (feature extraction) step that itself can be ill-posed and even further amplify errors. In order to reduce the error propagation of the first step, regularization strategies are usually applied. The choice of a suitable regularization strategy strongly depends on the particular situation and on the available prior information about the sought object . However, the recovery of requires different prior knowledge than feature extraction. This dismatch can lead to a substantial loss of performance in the feature detection step.

  5. In order to overcome the limitations mentioned in the remark above, image reconstruction and edge detection were combined in [louis2, louis1], where explicit formulas for estimating the edge map have been derived using the method of approximate inverse. This approach is also based on the dual convolution identity 3 and is closely related to the standard filtered backprojection (FBP) algorithm. However, this approach is not applicable to the case of undersampled data, since [louis1, louis2] employ the dual convolutional identity 3 and calculate the reconstruction filters of the form . In this calculation, in order to achieve a good approximation of the integral in 2.2, a properly sampled Radon data is required.

To overcome the limitations mentioned in the last remark above, we derive a novel framework for feature reconstruction in the next subsection (based on the forward convolutional identity 3) that does not make use of the continuous backprojection and, hence, can be applied to more general situations.

3.1 Proposed feature reconstruction

Our proposed framework for solving the feature reconstruction Problem 1 is based on the forward convolution identity 2 stated in Lemma 2.1. Because the convolution on the right-hand side of 2 acts only on the second variable, the convolution identity is not affected by the subsampling in the angular direction. Therefore, we have

(3.1)

Formally, the solution of (3.1) takes the form . If the data is properly sampled, this can be accurately and efficiently computed by applying the FBP algorithm to the filtered CT data . In this context, the data filter needs to be precomputed (from a given feature extraction filter ) in a filter design step. However, if the data is not properly sampled, the equations (3.1) are underdetermined and, in this case, FBP doesn’t produce accurate results, cf. [natterer, juergentodd]. In order to account for data incompleteness and stably approximate the feature map , a-priori information about the specific feature kernel or the feature map needs to be integrated into the reconstruction procedure. As a flexible way for doing this, we propose to approximate the inverse by the following variational regularization scheme:

(3.2)

Here denotes the noisy (semi-discrete data) and is a regularization (penalty) term.

Example 3.2.
  1. [itemsep=1ex]

  2. Image reconstruction: Here, the feature extraction filter is chosen as an approximation to the Delta distribution. For example, as with

    (3.3)

    being the Gaussian kernel. Another way of choosing for reconstruction purposes is through ideal low-pass filters

    , that are defined in the Frequency domain via

    , where , denotes a ball in with radius , and

    is the characteristic function of the set

    . It can be shown that in both cases as . These filters and its variants are often used in the context of the FBP algorithm.

  3. Gradient reconstruction: Here is chosen as a partial derivative of an approximation of the Delta distribution. For example, as with , . In this way, one obtains an approximation of the gradient of via

    where in the last equation above we applied the convolution componentwise. Such approximations of the gradient are for example used inside the well-known Canny edge detection algorithm [canny1986computational].

  4. Laplacian reconstruction: Analogously to the gradient approximation, is chosen to be the Laplacian of an approximation to the Delta distribution. A prominent example, is the Laplacian of Gaussian (LoG), i.e., , also known as the Marr-Hildreth-Operator. This operator is also used for edge detection, corner detection and blob detection, cf. [marr1980theory].

Depending on the problem at hand, there are several different ways of choosing the regularizer . Prominent examples in the case of image reconstruction include total variation (TV) or the -norm (possibly in the context of some basis of frame expansion). For the reconstruction of the derivatives (or edges in general), we will use the -norm as regularization term because derivatives of images can be assumed to be sparse and because the problem (3.2) can be efficiently solved in this case.

3.2 Filter design

The first step in our framework is filter design for (3.2). That is, given a feature extraction kernel , we first need to calculate the corresponding filter for the CT data, cf. (3.1). In our setting, filter design therefore amounts to the evaluation of the Radon transform of . In contrast to our approach, the filter design step in [louis2] consists in the calculation of a solution of the dual equation , given the feature extraction filter . As discussed above, the latter case requires full data and might be computationally more involved. From this perspective, filter design required by our approach offers more flexibility and can be considered somewhat simpler.

We now discuss some of the Examples 3.2 in more detail and calculate the associated CT data filters . In particular, we focus on the Gaussian approximations of the Delta distributions stated in (3.3). In a first step we compute the Radon transform of a Gaussian.

Lemma 3.3.

The Radon transform of the Gaussian , defined by (3.3), is given by

(3.4)

Since the Gaussian converges to the Delta distribution as , the smoothed version constitutes an approximation to for small values of . In order obtain approximations to partial derivatives of , we note that . Hence, using the feature extraction filters , the Problem 1 amounts to reconstructing partial derivatives of . Using this observation together with Lemma 3.3 and the property 4, we can explicitly calculate data filters used in different edge reconstruction algorithms (such as Canny or for the Marr-Hildreth-Operator).

Proposition 3.4.

Let the Gaussian be defined by (3.3).

  1. [itemsep=1ex]

  2. Gradient reconstruction: For the feature extraction filter the corresponding data filter is given by

    (3.5)

    Note that in (3.5), the notation refers to a vector valued function that is defined by a componentwise application of the Radon transform (cf. Example 3.2, No. 2).

  3. Laplacian reconstruction: For the feature extraction filter the corresponding data filter is given by

    (3.6)

From Proposition 3.4 we immediately obtain an explicit reconstruction formula for the approximate computation of the gradient and of the Laplacian of :

Both of the above formulas are of FBP type and can be implemented using the standard implementations of the FBP algorithm with a modified filter. This approach has the advantage that only one data filtering step has to performed, followed by the standard backprojection operation.

In order to derive FBP-filters for the gradient and Laplacian reconstruction, let us first note that , where the operator acts on the second variable and is defined in the Fourier domain by for , cf. [natterer]. Now, using the relations for the Fourier transform in 1D, , and , together with

(3.7)

we obtain the following result.

Proposition 3.5.

Let the FBP-filters and be given in the Fourier domain (componentwise) by

(3.8)

and

(3.9)

where and . Then, for , we have

(3.10)

Since the FBP algorithm is a regularized implementation of (cf. [natterer]), a standard toolbox implementation could be used in practice in order to compute and . To this end, one only needs to use the modified filters for the FBP, provided in (3.8) and (3.9), instead of the standard FBP filter (such as Ram-Lak). Again, let us emphasize that the reconstruction formulae (3.10) can only be used in the case of properly sampled CT data. If the CT data does not satisfy the sampling requirements, e.g., in case of angular undersampling, this FBP algorithm will produce artifacts which can substantially degrade the performance of edge detection. In such cases, our framework (3.2) should be used in combination with a suitable regularization term. In the context of edge reconstruction, we propose to use -regularization in combination with -regularization. This approach will be discussed in the next section.

So far, we constructed data filters for the approximation of the gradient and Laplacian in the spatial, cf. Proposition 3.4, and derived according FBP filters in the Fourier domain in Proposition 3.5. In a similar fashion, one can derive various related examples by replacing the Gaussian by feature kernels whose Radon transform is known analytically. Another way of obtaining practically relevant data filters (for a wide class of feature filter) is to derive expressions for the data filters in the Fourier domain (i.e., filter design in the Fourier domain). In the following, we provide two basic examples for filter design in the Fourier domain. To this end, we will employ the Fourier slice theorem, cf. Lemma 2.1, 1.

Remark 3.6.
  1. [itemsep=1ex]

  2. Lowpass Lalpacian: The Laplacian of the ideal lowpass is defined as

    where is the bandwidth of . Using the property 5, we get . By the Fourier slice theorem, we obtain

    Hence, the associated data filter is given by

    (3.11)

    Because is -band-limited, the convolution with the filter (3.11) can be discretized systematically whenever the underlying image is essentially band-limited. To that end, assume that the function has bandwidth . Then has bandwidth as well (with respect to the second variable), and therefore the continuous convolution can be exactly computed via discrete convolution. Using discretization (2.3) and taking we obtain from (3.11) the discrete filter

    (3.12)

    According to one-dimensional Shannon sampling theory, we compute via discrete convolution with the filter coefficients given in (3.12).

  3. Ram-Lak type filter: Consider the feature extraction filter

    where . Note that for , we have , since in this case . Hence, we consider the case . In a similar fashion as above, we obtain

    (3.13)

    Evaluating at , we get

    (3.14)

    Again, we can evaluate via discrete convolution with the filter coefficients (3.14).

Finally, let us note that there are several other examples for feature reconstruction filters for which one can derive explicit formulae of corresponding data filters in a similar way as we did in this section. For example, in the case of approximation of Gaussian derivatives of higher order or for band-limited versions of derivatives.

4 Numerical results

In our the numerical experiments, we focus on the reconstruction of edge maps. To that end, we use our framework (3.2) in combination with feature extraction filters that we have derived in Proposition 3.4 and in Remark 3.6. Since the gradient and the Laplacian of an image have relatively large values only around edges and small values elsewhere, we aim at exploiting this sparsity and, hence, use a linear combination as regularizer in (3.2). The resulting minimization problem then reads

(4.1)

If , this approach reduces to the regularization which is known to favor sparse solution. If , the additional -term increases smoothness of the recovered edges. In order to numerically minimize (4.1), we use the fast iterative shrinkage-thresholding algorithm (FISTA) of [fista]. Here, we apply the forward step to and the backward step to . The discrete norms are defined by and the discrete Radon transform

is computed via the composite trapezoidal rule and bilinear interpolation. The adjoint Radon transform

is implemented as a discrete backprojection following [natterer].

4.1 Reconstruction of the Laplacian feature map

We first investigate the feasibility of the proposed approach for recovering the Laplacian of the initial image. For our first experiment, we use a phantom image, which is defined as a characteristic function of the union of three discs and has the size with , cf. Figure 0(a). Since, according to the sampling condition (2.5), full aliasing free angular sampling requires samples in the -variable, we computed tomographic data at equally spaced signed distances and at equally spaced directions in . This data is properly sampled in the -variable, but undersampled in the angular variable , cf. Figure 0(b).

From this data, we computed the approximate Laplacian reconstruction, shown in Figure 0(c), using the standard FBP algorithm in combination with the LoG-filtered data that we computed in a preprocessing step using the LoG data filter from Proposition 3.4. It can be clearly observed that FBP introduces prominent undersampling artefacts (streaks), so that many edges in the calculated feature map are not related to the actual image features. This shows, that the edge maps computed by FBP (from undersampled data) can include unreliable information and even falsify the true edge information (since artefacts and actual edges superimpose). In a more realistic setup, this could be even worse, since artefacts may not be that clearly distinguishable from actual edges.

Figure 4.2 shows reconstructions of feature maps from noise-free CT data that we computed using our framework (4.1) for three different choices of feature extraction filters and for two different sets of regularization parameters. The first row of Figure 4.2 shows reconstructions with and using 1000 iterations of the FISTA algorithm, whereas the second row shows reconstructions that were computed using an additional -term with and using 500 iterations of the FISTA algorithm. In contrast to the FBP-LoG reconstruction (shown in Figure 0(c)), the undersampling artefacts have been removed in all cases. As expected, the use of -regularization without an additional -smoothing (shown in first row) produces sparser feature maps, as opposed to reconstruction shown in the second row. However, we also observed, that the iterative reconstruction based only on the -minimization (without the -term) sometimes has trouble reconstructing the object boundaries properly. In fact, we found that a proper reconstruction of boundaries is quite sensitive to the choice of the -regularization parameter. If this parameter is chosen too large, we observed that the boundaries could be incomplete or even disappear. Since the -regularization parameter controls the sparsity of the reconstructed feature map, this observation is actually not surprising. By including an additional -regularization term, the reconstruction results become less sensitive to the choice of regularization parameters.

[width=height=]Bilder/phantom

(a) Phantom

[width=height=]Bilder/data

(b) CT data

[width=height=]Bilder/del-fbp

(c) FBP-LoG reconstruction
Figure 4.1: Reconstruction of the Laplacian feature map using FBP. The phantom image of size consisting of a union of three discs (subfig:phantom) and the corresponding angularly undersampled CT data, measured at 40 equispaced angles in and properly sampled in the -variable with equispaced samples (subfig:data). Subfigure (subfig:fbp-log) shows the Laplacian of Gaussian (LoG) reconstruction using the standard FBP algorithm. It can be clearly observed that FBP introduces prominent streaking artefacts that are due to the angular undersampling.

[width=height=]Bilder/log

(a) LoG: ,

[width=height=]Bilder/lowpass

(b) Lowpass:

[width=height=]Bilder/ramlak

(c) Ram-Lak:

[width=height=]Bilder/logH1

(d) LoG:

[width=height=]Bilder/lowpassH1

(e) Lowpass:

[width=height=]Bilder/ramlakH1

(f) Ram-Lak:
Figure 4.2: Reconstruction of Laplacian feature maps using our framework. This figure shows reconstructions of feature maps from noise-free CT data that we computed using our framework (4.1) for three different choices of feature extraction filters and for two different sets of regularization parameters. Here, LoG refers to (3.6), Lowpass to (3.12), and Ram-Lak to (3.14). The first row shows reconstructions with and using 1000 iterations of the FISTA algorithm, whereas the second row shows reconstructions that were computed using an additional -term with and using 500 iterations of the FISTA algorithm. In contrast to the FBP-LoG reconstruction (shown in Figure 0(c)), the undersampling artefacts have been removed in all cases.

[width=height=]Bilder/reconstnoisenol1regul.png

(a) Ram-Lak:

[width=height=]Bilder/reconstdreiecknoise5000n40.png

(b) Ram-Lak:

[width=height=]Bilder/reconstdreiecknoiseH15000n40.png

(c) Ram-Lak:
Figure 4.3: Reconstructions of Laplacian feature maps from noisy data. Reconstruction in (subfig:noisy ramlak1) were calculated using only -regularization, in (subfig:noisy ramlak2) using only -regularization, and in (subfig:noisy ramlak3) using combined and -regularization.

In order to simulate the real world measurements more realistically, we added noise to the CT data that we used in the previous experiment. Using this noisy data, we calculated reconstructions via (4.1) in combination with the Ram-Lak type filter (3.14) using three different sets of regularization parameters and 1000 iterations of the FISTA algorithm in each case. The reconstruction using the parameters and (i.e., only -regularization was applied) is shown in Figure 2(a). The reconstruction in Figure 2(b) uses only regularization, i.e., and , and the reconstruction in Figure 2(c) applies both regularization terms with . In both reconstructions shown in Figure 2(b) and 2(c), the recovered features are much more apparent than for pure -regularization. As in the noise-free situation, we observe that the (pure) -regularization might generate discontinuous boundaries, whereas the combined --regularization results in smoother and (seemingly) better represented edges.

4.2 Edge detection

One main application of our framework for the reconstruction of approximate image gradients or approximate Laplacian feature maps is in edge detection. Clearly, feature maps that contain less artifacts can be expected to provide more accurate edge maps.

For this experiment, we used a modified phantom image that is shown in Figure 3(a). In contrast to the previously used phantom, this image includes also weaker edges that are more challenging to detect. For this phantom, we generated CT data using the same sampling scheme as in our first experiment (Section 4.1) and computed the LoG-feature maps using the FBP approach (cf. Figure 3(b)) and using our approach (cf. Figure 3(c)) with , , and iterations of the FISTA algorithm for (4.1). Subsequently, we generated corresponding binary edge maps by extracting the zero crossings of these LoG-feature maps (cf. Figure 3(d) and 3(e)) by using the MATLABs edge functions. Note that this procedure is one of the standard edge detection algorithm, known as the LoG edge detector, cf. [marr1980theory]

. For both methods, we took standard deviation

for the application of the Gaussian smoothing and threshold for the detection of the zero crossings. As can be clearly seen from the results, the edge detection based on our approach (cf. Figure 3(d)) is able to detect also the weaker edges inside the large disc. In contrast, edge detection in combination with FBP-LoG featuremap was not able to detect the edge set correctly due to strong undersampling artefacts.

[width=height=]edge/01_phantom.png

(a) Modified phantom

[width=height=]edge/03_LOG_FBP_phantom

(b) FBP-LoG

[width=height=]edge/06_Thresh_LOG_FBP_phantom

(d) FBP-LoG: edge map

[width=height=]edge/05_Thresh_L1_recon_phantom

(e) LoG: edge map
(c) LoG:

[width=height=]edge/02_FBP_phantom

Figure 4.4: LoG edge detection. The modified phantom image (subfig:mod phantom) includes also weaker edges that are more challenging to detect. Subfigures (3(b)) and (3(c)) show reconstructions of the LoG feature maps, that were generated using the FBP algorithm and our approach, respectively. The corresponding binary edge masks generated by the LoG edge detector are shown in (3(d)) and (3(e)).

[width=height=]lotus/lotus_fbp.png

(a) FBP reconstruction

[width=height=]lotus/lotus_sigma6_Gmag_fbp

(b) FBP-grad

[width=height=]lotus/lotus_sigma6_Gmag_l1reg

(c) grad: ,

[width=height=]lotus/lotus_sino.png

(d) CT data of a lotus

[width=height=]lotus/lotus_sigma6_Emap_fbp

(e) FBP-grad: edge map

[width=height=]lotus/lotus_sigma6_Emap_l1reg

(f) grad: edge map
Figure 4.5: Canny edge detection from the lotus data set. Rebinned CT data of a lotus root (subfig:lotus sino) (cf. [bubba2016tomographic]) and the corresponding FBP reconstruction (subfig:lotus fbp) from evenly distributed angles in . Magnitude of the smooth gradient map computed using the FBP algorithm (subfig:lotus fbp gradient) and using our approach (subfig:lotus log gradient). The corresponding edge detection results using the Canny algorithm are shown in (subfig:lotus fbp edge) and (subfig:lotus log edge), respectively.

In our last experiment we present edge detection results for real CT data of a lotus [bubba2016tomographic]. Note that similar reconstructions were presented in [frikel2021combining]. The lotus data have been rebinned and downsampled to signed distances and directions. The Gaussian gradient feature map was computed in two ways: First, by applying FBP to filtered CT data with the data filter (3.5), cf. Figure 4(b). Second, by using our approach (3.2) with and , and by applying 50 iterations of the FISTA algorithm, cf. Figure 4(c). The resulting image size is . The standard deviation for the Gaussian smoothing was chosen as and for the Canny edge detection we used the same lower threshold and upper threshold . In order to calculate binary edge maps (shown in Figures 4(e) and 4(f)), we used the Canny edge detector (cf. [canny1986computational]) in combination with the pointwise magnitude of the Gaussian gradient maps . Again, it was observed that the calculation of the Gaussian gradient map using our approach leads to more reliable edge detection results.

5 Conclusion

In this paper, we proposed a framework for the reconstruction of features maps directly from incomplete tomographic data, without the need of reconstructing the tomographic image first. Here, a features map refers to the convolution where is a given convolution kernel and is the underlying object. Starting from the forward convolution identity for the Radon transform, we introduced a variational model for feature reconstruction, that is formulated using the discrepancy term and a general regularizer . In contrast to existing approaches, such as [louis2, louis1], our framework does not require full data and, due to the variational formulation, also offers a flexible way for integrating a priori information about the feature map into the reconstruction. In several numerical experiments, we have illustrated that our method can outperform classical feature reconstruction schemes, especially, if the CT data is incomplete. Although we mostly focused on reconstruction of feature maps that are used for edge detetion purposes, our framework can be adapted for a wide range of problems. A rigorous convergence analysis of the presented scheme remains an open issue. Another directions of further research may may include the extension of the proposed approach to non-sparse, non-convolutional features and generalization to other types of tomography problems. Also, multiple feature reconstruction (similar to the method [jointmarkus, zangerl2020multi]) seems to be an interesting future research direction.

Acknowledgment

The work of the M.H was supported by Austrian Science Fund (FWF) project P 30747-N32. The contribution by S.G is part of a project that has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 847476. The views and opinions expressed herein do not necessarily reflect those of the European Commission.

References