Entropy Regularized Iterative Weighted Shrinkage-Thresholding Algorithm (ERIWSTA): An Application to CT Image Restoration

1 Introduction

For sparse compressed signals, Donoho $e t a l .$ [1] proposed a compressive sensing theory that enables efficient data sampling at a much lower rate than the requirements, which can be modeled as follows in its standard formulation.

$N o t a t i o n s$ : In this paper, the matrices are represented in capital letters. For a matrix $A$ , $A_{* i}$ , $A_{i *}$ and $A_{i j}$ denote the $i - t h$ column, the $i - t h$ row and $(i, j) - t h$ element of $A$ , respectively; the $∥ {˙ ∥}_{i}$ represents the $i$

-norm of a vector. All the vectors are column vectors unless transposed to a row vector by a prime superscript

T

Compressive sensing can be formulated as:

b = A x + ϵ

(1)

where $x \in R^{n}$ is an unknown vector, $b \in R^{m}$ is an observed vector, $A \in R^{m \times n}$ is called the compressive sensing matrix (usually $m << n$ ), and $ϵ$ is the unknown disturbance term or noise. Obviously, this is an underdetermined system of equations that does not have a sole solution. The least-squares method is usually used to solve the problem.

min x \frac{1}{2} ∥ A x - b ∥_{2}^{2}

(2)

To suppress overfitting, some scholars [2, 3, 4, 5] added the $L_{0}$ -norm regularizer to introduce sparse prior information.

min x \frac{1}{2} ∥ A x - b ∥_{2}^{2} + β ∥ x ∥_{0}

(3)

where $∥ x ∥_{0}$ denotes the number of nonzero components of $x$ and $β > 0$

is a hyperparameter to control the tradeoff between accuracy and sparsity. Many methods have been developed to solve this problem, such as the penalty decomposition method

[6], iterative hard threshold method [7], fixed-point continuation method (FPC) [9], approximate gradient homotopy method (PGH) [10] and reweighted

L_{0}

minimization method [11, 12].

However, Eq. (3) is an NP-hard optimization problem [13], which is highly discrete so that it is difficult to solve using a precise algorithm. Thus, we need to seek an effective approximation solution for this problem. The $L_{1}$ -norm regularizer is introduced as a substitute for that of the $L_{0}$ -norm. Such an approximation can be traced back to a wide range of fields, such as seismic traces [14], sparse signal recovery [15], sparse model selection in statistics (LASSO) [16], and image processing [17]. Many scholars have attempted to find the optimal solution to the following problem:

min x \frac{1}{2} ∥ A x - b ∥_{2}^{2} + β ∥ x ∥_{1}

(4)

It is a convex continuous optimization problem with a sole nondifferentiable point ( $x = 0$ ), which can usually be transformed into a second-order cone programming problem and then solved by methods such as interior-point methods. However, in large-scale problems, due to the high algorithmic complexity, the interior-point method is very time-consuming. Based on this, many researchers have solved the problem through simple gradient-based methods. Among them, the iterative shrinkage-thresholding algorithm (ISTA) proposed by Chambolle $e t a l .$ [18, 19] has attracted much attention. ISTA updates $x$ through a shrinkage/soft threshold operation in each iteration.

x^{k + 1} = s o f t_{β t} [x^{k} - 2 t A^{T} (A x^{k} - b)]

(5)

where $k$ represents the $k$ -th iteration, $t$ is an appropriate stepsize and $s o f t$ is the soft threshold operation function.

s o f t_{θ} (x_{i}) = s i g n (x_{i}) (∥ x_{i} ∥ - θ)

(6)

Recently, the iteratively weighted shrinkage-thresholding algorithm (IWSTA) has attracted much interest compared with ISTA, which outperforms their unweighted counterparts in most cases. In these methods, decision variables and weights are optimized alternatingly, or decision variables are optimized under heuristically chosen weights. It can be written as:

min x, w \geq 0 \frac{1}{2} ∥ A x - b ∥_{2}^{2} + β n \sum i = 1 w_{i} | x_{i} |_{1}

(7)

The method assigns different weights to each component of $x$ in the iterative process and then updates $x$ . In this way, each subproblem is convex and easy to solve. Many algorithms have been developed to solve it. For example, the iterative support detection (ISD) method [20] assigns a weight of 0 to components in the support set and a weight of 1 to the other components during iteration, in which the support set at each iteration consists of all components whose absolute value is greater than the threshold. Zhao $e t a l .$ [12] proposed a new method to calculate the optimal $w$

by the duality of linear programming based on the property of weighted range space. It alternately solves the weighted original problem with fixed weights to obtain a new solution

x

, and then it solves the duality problem to obtain a new weight

w

. More variants are available in [11, 21] and its references. The details of some examples are listed in Tab. 1.

Author	Termed	Weights	Min.	Max.	Regularizer
Chambolle $e t a l .$ [18]	ISTA	1	1	1	$n \sum i = 1 \| x_{i} \|$
Candes $e t a l .$ [22]	IRL1	$\frac{1}{\| x_{i}^{k - 1} \| + δ}$	0	$\frac{1}{δ}$	$n \sum i = 1 log (\| x_{i} \| + δ)$
Foucart $e t a l .$ [23]	WLP	$\frac{1}{(\| x_{i}^{k - 1} \| + δ)^{1 - p}}$	0	$\frac{1}{δ^{1 - p}}$	$\frac{1}{p} n \sum i = 1 (\| x_{i} \| + δ)^{p}$
Wipf $e t a l .$ [24]	NW4	${\frac{1 + (\| x^{k - 1} \| + δ)^{p + 1}}{(\| x^{k - 1} \| + δ)}}^{p + 1}$	0	1	$n \sum i = 1 (\| x_{i} \| - \frac{1}{(x_{i} + δ)^{p}})$

Table 1: Variants of weighted method.

There is a drawback for the above methods: the weights do not meet the usual definition of weights, and their sum is one, which leads them to be distributed in a very large range (see Tab. 1). Such weights are difficult to explain and can lead to an inaccurate result.

This paper proposes a new IWSTA type, called entropy regularized IWSTA (ERIWSTA), which obtains easily computable and interpretable weights. The weights automatically fall in the range of [0, 1], and the summation is one so that they can be considered a probability of the contribution of each attribute to the model. This is achieved by adding an entropy regularizer to the cost function and then using the Lagrange multiplier method to solve the problem. Experiments are executed for CT image restoration, and the results show that the proposed algorithm performs better in terms of both convergence speed and restoration accuracy compared with some state-of-the-art methods.

2 Methodology

The main idea of the IWSTA type algorithms is to define a weight for each attribute based on the current iteration $x^{k}$ and then use them to obtain a new $x$ . In this section, we introduce an entropy regularizer to the cost function and obtain the following optimization model:

$min$	$Φ_{β, γ} (x, w) = F (x) + β G_{γ} (x, w)$
$s . t .$	$w_{i} \geq 0, n \sum i = 1 w_{i} = 1$
$w h e r e$	$F (X) = \frac{1}{2} ∥ A x - b ∥_{2}^{2}$	(8)
	$G_{γ} (x, w) = n \sum i = 1 w_{i} \| x_{i} \| + γ n \sum i = 1 w_{i} ln w_{i}$	(8)

where $γ \geq 0$ is a given hyperparameter.

As can be seen, while we do not use the entropy regularizer, $w$ can easily be solved as $w_{i} = 1$ if $x_{i} = arg min {| x_{1} |, . . ., | x_{n} |}$ , or 0 otherwise¹¹1The update rule can be easily explained by an example as

min {4, - 1, 5} = m i n

4 w_{1} - 1 w_{2} + 5 w_{3}

s . t .

w_{1}, w_{2}, w_{3} \geq 0

w_{1} + w_{2} + w_{3} = 1

The solution is

w_{1} = 0

w_{2} = 1

and

w_{3} = 0

, in which

w_{2}

corresponds to the minimum value of {4, -1, 5}. It is very similar to the computation of the weights in the k-means algorithm.

. It shows a simple fact that only one element of

w

is 1, and the others are 0, which is grossly incompatible with the actual problem. Then, we add the negative entropy of the weights to measure the uncertainty of weights and stimulate more attributes to help signal reconstruction because it is well known that

\sum_{i = 1}^{n} w_{i} ln w_{i}

is minimized in information theory when

w_{1} = w_{2} = . . . = w_{n}

(9)

As follows, we will alternatively solve $w$ and $x$ in Eq. (2).

2.1 Update rule for $w$

To solve $w$ , we introduce the Lagrange multipliers $λ$ and then obtain the following Lagrange function. Note that $F (x)$ is a constant with respect to $W$ , so we only construct the Lagrange function on $G (x)$ .

L (w, λ) = G_{γ} (x, w) + λ (n \sum i = 1 w_{i} - 1),

(10)

Set the partial derivative of $L (w, λ)$ with respect to $w_{i}$ and $λ$ to zero and then obtain the following two equations.

	$\frac{\partial L (w, λ)}{\partial w_{i}}$	$=$	$\| x_{i} \| + γ (1 + ln w_{i}) + λ = 0$		(11)
	$\frac{\partial L (w, λ)}{\partial λ}$	$=$	$n \sum i = 1 w_{i} - 1 = 0$		(12)

From Eq. (11), we know that

w_{i} = exp (- \frac{λ}{γ}) exp (- \frac{| x_{i} |}{γ})

(13)

Substituting Eq. (13) into Eq. (12), we have

n \sum i = 1 w_{i} = exp (- \frac{λ}{γ}) n \sum i = 1 exp (- \frac{| x_{i} |}{γ}) = 1

(14)

It follows that

exp (- \frac{λ}{γ}) = \frac{1}{\sum_{i = 1}^{n} exp (- \frac{| x_{i} |}{γ})}

(15)

Substituting this expression to Eq. (13), we obtain that

w_{i} = \frac{exp (- \frac{| x_{i} |}{γ})}{\sum_{l = 1}^{n} exp (- \frac{| x_{l} |}{γ})}

(16)

Such weights certainly meet the constraints that $w_{i} \geq 0$ and $\sum_{i = 1}^{n} w_{i} = 1$ .

2.2 Update rule for $x$

Inspired by the work of ISTA [8], a similar approach was adopted for the iterative update of $x$ . The construction of a majorization is an important step in obtaining the updating rule.

Definition 2.1

(Majorization) Denote $ψ (x | x^{k})$ as a majorization for $Ψ (x)$ at $x^{k}$ (fixed) if $ψ (x^{k} | x^{k}) = Ψ (x^{k})$ and $ψ (x | x^{k}) \geq Ψ (x)$ .

Clearly, $Ψ (x)$ is nonincreasing under the updating rule $x^{k + 1} = {min}_{x} ψ (x | x^{k})$ because

Ψ (x^{k + 1}) \leq ψ (x^{k + 1} | x^{k}) \leq ψ (x^{k} | x^{k}) = Ψ (x^{k})

(17)

Then, we can construct the majorization for $F (x)$ .

Proposition 2.1

Obviously, $F (x)$ is a Lipschitz continuous and differentiable convex function, which has a majorization function at fixed current iteration $x^{k}$ as

f (x, x^{k}) = F (x^{k}) + [\nabla F (x^{k})]^{T} (x - x^{k}) + \frac{L}{2} ∥ x - x^{k} ∥_{2}^{2}

(18)

where $L$

is larger than or equal to the maximum eigenvalue of

A^{T} A

Proof 2.1

It is well-known that

F (x) = \frac{1}{2} ∥ A x - b ∥_{2}^{2} = F (x^{k}) + [\nabla F (x^{k})]^{T} (x - x^{k}) + \frac{}{1} 2 (x - x^{k})^{T} A^{T} A (x - x^{k})

(19)

We compare $F (x)$ and $f (x, x^{k})$

and find that only the last terms are different. By singular value decomposition (SVD) of a symmetric definite matrix, we know that

A^{T} A = Q^{T} Σ Q

, in which

Q

is an orthogonal matrix consisting of all eigenvectors and

Σ

is diagonal consisting of all eigenvalues. Let

z = x - x^{k}

, then

z^{T} (A^{T} A) z = z^{T} Q^{T} Σ Q z \leq L ∥ Q z ∥_{2}^{2} = L ∥ z ∥_{2}^{2}

(20)

And it is also certain that $z^{T} A^{T} A z = L ∥ z ∥_{2}^{2} = 0$ if $x = x^{k}$ . Thus, the proof is established.

Now, we obtain the majorization for the cost function $Φ (x, w)$ on $x$ .

ϕ (x, x^{k}) = f (x, x^{k}) + β G_{γ} (x, w)

(21)

which can be reorganized as

	$ϕ (x, x^{k})$	$=$	$\frac{L}{2} ∥ x - [x^{k} - \frac{1}{L} \nabla F (x^{k})] ∥_{2}^{2} + β G_{γ} (x, w)$		(22)
		$=$	$n \sum i = 1 {\frac{L}{2} ∥ x_{i} - [x^{k} - \frac{1}{L} \nabla F (x^{k})]_{i} ∥_{2}^{2} + β w_{i} \| x_{i} \|} + c o n s t a n t$		(22)

We find that the variables of the majorization are separable such that their minimizations can be easily obtained on each $x_{i}$ , respectively, as follows:

x_{i}^{k + 1} = s o f t_{β t w_{i}} [x^{k} - 2 t A^{T} (A x^{k} - y)]

(23)

Figure 1: The original and noisy head phantom images. $(a)$ head phantom with 256×256 pixels; $(b)$ and $(c)$ blurred image with a 5 $\times$ 5 uniform kernel and additive Gaussian noise with $σ = 10^{- 2}$ and $σ = 10^{- 3}$ .

Figure 2: 3D profile of $β$ and $λ$ on MAE with different Gaussian noise levels: (a) $σ = 10^{- 2}$ and (b) $σ = 10^{- 3}$ .

Termed	$β$	$γ$	$δ$	MAE
ISTA	$10^{- 3}$	$-$	$-$	$5.312077 * 10^{- 7}$
WLP	$10^{- 5}$	$10^{- 10}$	$10^{- 3}$	$5.228672 * 10^{- 7}$
NW4	$10^{- 5}$	$10^{- 2}$	$10^{- 3}$	$5.410231 * 10^{- 7}$
IRL1	$10^{- 5}$	$-$	$10^{- 3}$	$5.228672 * 10^{- 7}$
ERIWSTA	$10^{2}$	$10^{- 2}$	$-$	$5.218246 * 10^{- 7}$

Table 2: The optimal MAE value and corresponding hyperparameter (Gaussian noise with

σ = 10^{- 2}

Termed	$β$	$γ$	$δ$	MAE
ISTA	$10^{- 3}$	$-$	$-$	$5.122013 * 10^{- 7}$
WLP	$10^{- 5}$	$10^{- 5}$	$10^{- 3}$	$5.018339 * 10^{- 7}$
NW4	$10^{- 5}$	$10^{- 2}$	$10^{- 3}$	$5.410231 * 10^{- 7}$
IRL1	$10^{- 5}$	$-$	$10^{- 3}$	$5.018340 * 10^{- 7}$
ERIWSTA	$10^{2}$	$10^{- 2}$	$-$	$5.005524 * 10^{- 7}$

Table 3: The optimal MAE value and corresponding hyperparameter (Gaussian noise with

σ = 10^{- 3}

Figure 3: Cost function versus iteration number for different Gaussian noise levels: (a) $σ = 10^{- 2}$ and (b) $σ = 10^{- 3}$ .

Figure 4: MAE versus iteration number for different Gaussian noise levels: (a) $σ = 10^{- 2}$ and (b) $σ = 10^{- 3}$ .

Figure 5: After 30 iterations, the denoising results of ISTA, WLP, NW1, IRL1 and ERIWSTA with Gaussian noise with $σ = 10^{- 2}$ .

Figure 6: After 30 iterations, the denoising results of ISTA, WLP, NW1, IRL1 and ERIWSTA with Gaussian noise with $σ = 10^{- 3}$ .

Figure 7: Horizontal central profiles of the restored images with different Gaussian noise levels: (a) $σ = 10^{- 2}$ and (b) $σ = 10^{- 3}$ ..

Figure 8: Vertical central profiles of the restored images with different Gaussian noise levels: (a) $σ = 10^{- 2}$ and (b) $σ = 10^{- 3}$ .

3 Numerical experiments

Numerical experiments are provided to evaluate the performance of the proposed ERWISTA compared with ISTA, WLP, NW4 and IRL1 on the denoising problem of computed tomography (CT) images. All experiments are performed on an HP computer with a 2.5 GHz Intel(R) Core(TM) i7-4710MQ CPU with 12 GB of memory using MATLAB R2019a for coding. A simulated Shepp-Logan phantom with $256 \times 256$ pixels was used to evaluate the algorithm performance, which is usually used in CT image analysis. There are many advantages to using simulated phantoms, including prior knowledge of the pixel values and the ability to control noise. We blurred the image by using a uniform $5 \times 5$ kernel (applied by the MATLAB function "fspecial") and then added Gaussian noise by the following formula. We select $σ = 10^{- 2}$ and $10^{- 3}$ as examples of high and low noise levels for the following experiments.

x^{n o i s e} = x^{t r u e} + N (0, σ)

(24)

Fig. 1 shows the original and blurred-and-noisy images. Based on the good time-frequency localization characteristics of the wavelet transform, it can effectively distinguish high-frequency noise from low-frequency information. Therefore, the wavelet transform is used to reduce noise. The introduction of the wavelet matrix can also ensure the sparsity of the whole optimization algorithm. Without losing generality, let $A = P W$ , where $P$ is the predetermined system matrix indicating the blurring information and $W$ represents the second-order Haar wavelet matrix.

Mean absolute error (MAE) was used to measure the similarity to the true image. The value of the MAE was calculated by taking the average of the squared differences between the restored pixel values and the true pixel values.

M A E = \frac{1}{N} | | x^{r e s t o r a t i o n} - x^{t r u e} | |_{1}

(25)

3.1 Hyperparameter selection

To select the penalty hyperparameter $β$ and the entropy weighted hyperparameter $γ$ , we compare the MAE value after 100 iterations with respect to them from $10^{- 10}$ to $10^{10}$ . The results are shown in Fig. 2, demonstrating that ERIWSTA can achieve a consistently low MAE value over a wide range of $β$ and $γ$ , which displays its robustness.

We also quantitatively display the optimal MAE and corresponding hyperparameters of the algorithms in Tabs. 2 and 3. An interesting observation is that, regardless of low or high noise levels, the restoration accuracy of our algorithm is always better than the others. These optimal hyperparameters are also used in the following experiments.

3.2 Algorithmic performance

Fig. 3 displays the cost function of the algorithms. As can be seen, the proposed algorithm always have the fast convergence speed, which arrive at the stable status early.

Fig. 4 shows the MAE curves of the algorithms with respect to the number of iterations. The proposed ERIWSTA always has superior performance to the other algorithms, rapidly obtaining the minimum MAE value.

Figs. 5 and 6 indicate the denoising results with the given noise level. As can be seen, all of the algorithms achieve a similar image. Howver, Figs. 7 and 8 quantitatively compares the vertical profiles of the restored images with that of the true phantom in the central row and column. We can see that ERIWSTA follows the outline of the phantom more accurately than the other algorithms.

4 Conclusions

In this paper, a new IWSTA type, called ERIWSTA, is proposed to solve the linear inverse problem. An entropy weighted term is introduced to measure the certainty of the weights, and then the Lagrange multiplier method is used to obtain a simple solution. The experimental results on image restoration of a synthetic CT head phantom show that ERIWSTA can achieve convergence faster with fewer iterations and better restoration accuracy than the other algorithms.

However, as with many existing algorithms, our algorithm also involves two main hyperparameters ( $β$ and $γ$ ). In the future, we will focus on designing an automatic method to adjust these hyperparameters.

Acknowledgments

This work was supported by the Fundamental Research Funds for the Central Universities (N2019006 and N180719020).

References

[1] E. Candes, J. Romberg, T. Tao, Stable signal recovery from incomplete and inaccurate measurements, Commun. Pure Appl. Math. 59 (8) (2005) 1207-1223.
[2] M. Elad, M. A. Figueiredo, Y. Ma, On the role of sparse and redundant representations in image processing, Proc. IEEE 1998, 6 (2010) 972-982.
[3]
J. Wright, M. Yi, J. Mairal, G. Sapiro, T. S. Huang, S. Yan, Sparse representation for computer vision and pattern recognition, Proc. IEEE 1998, 6 (2010) 1031-1044.
[4] M. H. Alsuwaiyel, Sparse Coding and its Applications in Computer Vision, World Scientific, 1999.
[5]
M. Gong, J. Liu, H. Li, Q. Cai, L. Su, A multiobjective sparse feature learning model for deep neural networks, IEEE Trans. Neural Netw. Learn. Syst. 26 (12) (2015) 3263-3277.
[6] Z. Lu, Y. Zhang, Sparse approximation via penalty decomposition methods, SIAM J. Optimiz. 23 (4) (2012) 2448-2478.
[7] T. Blumensath, M. E. Davies, Iterative thresholding for sparse approximations, J. Fourier Anal. Appl. 14 (5) (2008) 629-654.
[8] A. Beck, M. Teboulle, A fast iterative shrinkage-thresholding algorithm for linear inverse problems, SIAM J. Imag. Sci 2 (1) (2009) 183-202.
[9] E. T. Hale, W. Yin, Y. Zhang, Fixed-point continuation for $L_{1}$ -minimization: Methodology and convergence, SIAM J. Optimiz. 19 (3) (2008) 1107-1130.
[10] L. Xiao, T. Zhang, A proximal-gradient homotopy method for the sparse least-squares problem, SIAM J. Optimiz. 23 (2) (2012) 1062-1091.
[11] Z. Xie, J. Hu, Rewighted $L_{1}$ -minimization for sparse solutions to underdetermined linear systems, 2013 6th Int. Congress Image Sig. Process. (CISP) 22 (2012) 1065–1088.
[12] M. K. Y.B.Zhao, A new computational method for the sparsest solutions to systems of linear equations, SIAM J. Optimiz. 25 (2) (2015) 1110-1134.
[13] B. K. Natarajan, Sparse approximate solutions to linear systems, SIAM J. Comput. 24 (2) (1995) 227–234.
[14] H. L. Taylor, S. C. Banks, J. F. Mccoy, Deconvolution with the $L_{1}$ -norm, Geophysics 44 (1) (1979) 39.
[15] S. S. Chen, D. Saunders, Atomic decomposition by basis pursuit, SIAM Rev. 43 (1) (2001) 129-159.
[16] R. J. Tibshirani, Regression shrinkage and selection via the lasso, J. R. Stat. Soc. B 58 (1) (1996) 267-288.
[17] P. V. Blomgren, J. D. Villasenor, B. Engquist, S. Osher, V. B. Bo, Total variation methods for restoration of vector valued images, SIAM J. Optimiz. 25 (2) 1110-1134.
[18] A. Chambolle, R. De Vore, N.-Y. Lee, B. Lucier, Nonlinear wavelet image processing: Variational problems, compression, and noise removal through wavelet shrinkage, IEEE Trans. Image Process. 7 (3) (1998) 319-335.
[19] M. Figueiredo, R. D. Nowak, An EM algorthim for wavelet-besed image restoration, IEEE Trans. Image Process. 12 (8) (2003) 906-916.
[20] Y. Wang, W. Yin, Sparse signal reconstruction via iterative support detection, SIAM J. Imag. Sci. 3 (3) (2009) 462-491.
[21] D. L., Donoho, For most large underdetermined systems of linear equations the minimal 1-norm solution is also the sparsest solution, Commun. Pure Appl. Math. 59 (6) (2006) 797-829.
[22] E. Candès, M. B. Wakin, S. P. Boyd, Enhancing sparsity by reweighted $L_{1}$ minimization, J. Fourier. Anal. Appl. 14 (5-6) (2008) 877-905.
[23] S.Foucart, M.Lai, Sparsest solutions of undetermined linear systems via $L_{q}$ -minimization for $0 < q <= 1$ , IEEE Trans. Image Process. 26 (3) (2009) 395-407.
[24] D. Wipf, S. Nagarajan, Iterative reweighted $L_{1}$ and $L_{2}$ methods for finding sparse solutions, IEEE J. Sel. Topics Signal Process. 4 (2) (2010) 317-329.

Entropy Regularized Iterative Weighted Shrinkage-Thresholding Algorithm (ERIWSTA): An Application to CT Image Restoration

Authors

Learning optimal nonlinearities for iterative thresholding algorithms

An Iterative Shrinkage Approach to Total-Variation Image Restoration

Learned Block Iterative Shrinkage Thresholding Algorithm for Photothermal Super Resolution Imaging

Robust lEarned Shrinkage-Thresholding (REST): Robust unrolling for sparse recover

An Online Plug-and-Play Algorithm for Regularized Image Reconstruction

A Trained Regularization Approach Based on Born Iterative Method for Electromagnetic Imaging

Iterative Shrinkage Approach to Restoration of Optical Imagery

1 Introduction