Noise Reduction

, which can be assumed to be two dimensional without a loss of generality, and impulse response 
$$ \varphi $$
. Let 
$$ \bf D\subset {{\mathbb{R}^2}} $$
be squared area of the size 
$$ D\times D $$
, and let 
$$ \varphi $$
be a continuous function of fast decay whose Fourier transform 
$$ \hat{\varphi}({\omega_x},{\omega_y}) $$
is supported inside the circle 
$$ \omega_x^2+\omega_y^2=\Omega_0^2 $$
, i.e., 
$$ \hat{\varphi}({\omega_1},{\omega_2})=0 $$
if 
$$ \omega_1^2+\omega_2^2>\Omega_0^2 $$
” src=”/wp-content/uploads/2016/08/A270527_1_En_8_Chapter_IEq00089.gif”></SPAN>. The output <SPAN id=IEq000810 class=InlineEquation><IMG alt= of the imaging system is modeled as a convolution




$$ F=f*\varphi, $$

(8.1)
where 
$$ f:\mathbf{D}\to \mathbb {R} $$
is the input signal. It is supposed that 
$$ f $$
is integrable on 
$$ \mathbf{D} $$
. Then 
$$ F $$
can be represented by the Fourier series



$$ F(x,y)=\sum\limits_{{m,k\in \mathbb{Z}}} {{{\hat{F}}_{{m,k}}}{{\mathrm{ e}}^{{\mathrm{ i}\delta (xm+yk)}}},} $$

(8.2)
where 
$$ {{\hat{F}}_{m,k }}={D^{-2 }}\hat{f}(\delta m,\delta k)\hat{\varphi}(\delta m,\delta k) $$
and



$$ \delta =2\pi /D. $$

(8.3)


Since the impulse response 
$$ \varphi $$
is band limited with the bandwidth 
$$ {\Omega_0} $$
, the series (8.2) can be truncated:



$$ F(x,y)=\sum\limits_{{-N\leq m,k\leq N}} {{{\hat{F}}_{m,k }}{{\mathrm{ e}}^{{\mathrm{ i}(xm+yk)\delta }}}} $$

(8.4)
with



$$ N=\left\lceil {\frac{D}{\Delta }} \right\rceil $$

(8.5)
and



$$ \Delta =\frac{{2\pi }}{{{\Omega_0}}} $$

(8.6)

That is, the output 
$$ F $$
is a trigonometric polynomial and, consequently, can be uniquely recovered from its samples (see, e.g., [1]) measured over the rectangular grid 
$$ {{\{({x_m},{y_k})|{x_m}=m\Delta, {y_k}=k\Delta \}}_{{-N\leq m,k\leq N}}}. $$
Hence, under certain idealizing assumptions, the continuous image can be identified with its samples.

We refer to the matrix 
$$ {{\{{F_{m,k }}\}}_{{-N\leq m,k\leq N}}} $$
as a digital image of the input signal 
$$ f $$
. In the following, we call the number 
$$ N $$
the resolution limit of the imaging system, meaning that 
$$ N\delta $$
is the maximal resolution that can be achieved with this system. Clearly, for the case considered here, the value 
$$ N\delta $$
coincides with the Nyquist frequency of the imaging system.1

Finally, we point out that extending the matrix 
$$ \{{F_{j,n }}\} $$
periodically with the period 
$$ M\times M $$
, 
$$ M=2N+1 $$
, we can write



$$ {{\hat{F}}_{l,m }}=\frac{1}{{{M^2}}}\sum\limits_{j=0}^{M-1 } {\sum\limits_{k=0}^{M-1 } {{F_{j,k }}{{\mathrm{ e}}^{{-2\pi \mathrm{ i}(jl+km)/M}}}} }, $$

(8.7)
which is the standard notation for the discrete 2D Fourier transform.



8.2 Image Processing


In practice, we do not have to deal with the image 
$$ F $$
but with its corrupted version 
$$ I $$



$$ I=F+\eta, $$

(8.8)
where 
$$ \eta $$
is the random value referred to as system noise. For an X-ray imaging system, another substantial noise component, which is due to the Poisson statistics and scattering of the X-ray quanta, is referred to as quantum noise (see, e.g., Chap. 9 of [2]). In view of the fact that the quantum noise is an inherit feature of the input signal 
$$ f $$
, the representation (8.8) can be replaced by



$$ I=(g+q)*\varphi +\eta, $$

(8.9)
where 
$$ q $$
is the input quantum noise and 
$$ g $$
is the desired signal. One of the most important problems of the image processing is to recover 
$$ g $$
from the measured image 
$$ I $$
. It is clear, that this can be done to some degree of accuracy only. The best possible approximation of 
$$ g $$
in terms of the mean squared error is given by



$$ Ag(x,y)=\sum\limits_{{-N\leq k,l\leq N}} {{{\hat{g}}_{k,l }}{{\mathrm{ e}}^{{\mathrm{ i}(xk+yl)\delta }}}}, $$

(8.10)
where, at the right, there is the Fourier series of 
$$ g $$
truncated up to the resolution limit of the given imaging system. Leaving aside details, we mention that the linear problem of recovering unknowns 
$$ {{\hat{g}}_{k,l }} $$
from 
$$ I $$
is sometimes referred to as the Wiener filtering.2

As it follows from (8.9), the problem of recovering 
$$ Ag $$
from 
$$ I $$
can be looked at as the one consisting of two distinct tasks: the suppression of the noisy component 
$$ q*\varphi +\eta $$
(noise reduction) and the deconvolution of blurred image 
$$ g*\varphi $$
(de-blurring). In the following, we use this discrimination in two tasks and describe some approaches used for the reduction of noise in the image.


8.3 Linear Filtering


The general strategy of this approach is to separate in the spectral description of the medical image the frequencies corresponding mainly to signals and those mainly corresponding to noise. Afterwards it is tried to suppress the power of those spectral components which are related to noise. This is done via appropriate weighting of the frequency components of image 
$$ I $$
. The result of the weighting is the image 
$$ wI $$
, given by



$$ wI(x,y)={D^{-2 }}\sum\limits_{{-N\leq l,m\leq N}} {{{\hat{w}}_{l,m }}{{\hat{I}}_{l,m }}{{\mathrm{ e}}^{{\mathrm{ i}(xl+ym)\delta }}}}. $$

(8.11)

The set 
$$ {{\{{{\hat{w}}_{l,m }}\}}_{{-N\leq l,m\leq N}}} $$
is called transfer matrix. The elements of transfer matrix are samples of the Fourier transform of the filter 
$$ w $$
, in terms of which the right-hand side of (8.11) can be written as the 2D convolution



$$ wI=I*w. $$

(8.12)

In the following, we call function 
$$ w\in {L_2}({{\mathbb {R}}^2}) $$
smoothing function if it is continuous and



$$ \hat{w}(0)=\int_{{{{\mathbb {R}}^2}}} {w(x,y)\mathrm{ d}x\mathrm{ d}y\ne 0}. $$

(8.13)

Due to (8.12) the kernel 
$$ w $$
acts upon the image 
$$ I $$
as a local weighted averaging, which is uniform for the whole image and is independent from its current space context. Linear filtering is perfect if the power spectrum of the image 
$$ I $$
can be separated in two disjoint regions, the one of which contains the spectral component of the useful signal and, the other, of the noise. We illustrate this apparent property in Fig. 8.1, where in the middle, there is an image corrupted with additive noise, the spectral power of which is concentrated outside the circle 
$$ \omega_x^2+\omega_y^2={\Omega^2} $$
. On the contrary, the spectral power of the useful signal (left) is concentrated inside this circle. The right image is the convolution of the corrupted image with the ideal filter, i.e., the filter, the transfer function of which is

A270527_1_En_8_Fig1_HTML.jpg


Fig. 8.1
Linear filtering of the band-limited image. Left: the image the power spectrum of which is inside the circle B. Middle: the same image with added noise, the power spectrum of which is outside the circle B. Right: convolution of the middle image with the ideal filter given by (8.14)




$$ \begin{array}{clclclcl}{ \hat{w}(\boldsymbol{ \upomega} )=\left\{{ 1,\quad |\boldsymbol{ \upomega} |\leq \Omega, 0,\quad \mathrm{ otherwise}}\right.}\end{array}$$

(8.14)

Clearly, in this case, the original image and the denoised one are identical. In fact, if the power spectrum of the useful signal is concentrated inside a bounded region of the frequency domain and the signal-to-noise ratio is big enough in this region, then the linear filtering is still excellent independently from the spectral distribution of the noise. The corresponding example is illustrated in Fig. 8.2. On the left of this figure, there is the same image as in Fig. 8.1 with added white noise and, on the right, the result of the filtering with Gaussian filter 
$$ {g_{\varOmega }} $$
, the transfer function of which is

A270527_1_En_8_Fig2_HTML.jpg


Fig. 8.2
Linear filtering of the band-limited image. Left: image from Fig. 8.1 with added white noise. Right: the convolution with Gaussian kernel




$$ {{\hat{g}}_{\Omega }}(\boldsymbol{ \upomega} )={{\mathrm{ e}}^{{-\frac{1}{2}{{{\left( {\frac{{|\boldsymbol{ \upomega} |}}{\Omega }} \right)}}^2}}}}. $$

(8.15)

Thus, the following short conclusion can characterize the linear filtering: it works well if applied to regular signals, i.e., such signals the power spectrum of which is concentrated within a bounded region of low frequencies and it vanishes fast outside this region. The filtering operation in this case is reduced to the suppression of high-frequency components of the image. If the Fourier transform of an image decays slowly as 
$$ |\boldsymbol{ \upomega} |\to \infty $$
, the linear filtering is not efficient anymore. Since the spectral power of such images is high also in the high-frequency sub-band, suppressing within this sub-band affects the sharpness of edges and destroys fine details of the image. This results in overall blurring of the image (see Fig. 8.3).

A270527_1_En_8_Fig3_HTML.jpg


Fig. 8.3
Blurring effect of the linear filtering. Left: X-ray image of a human lung specimen. Right: convolution with a Gaussian kernel

There is a variety of filters and corresponding smoothing kernels which are available. However, all rely on the same principle and show similar behavior.


8.4 Adaptive Nonlinear Filtering


In contrast to methods based on linear filtering, many advanced noise reduction techniques assume a control over the local amount of smoothing. Appropriately choosing a smoothing kernel at the current positions, it is possible to avoid unnecessary blurring in regions of high variation of the useful signal. Usually in practice the desired kernel is chosen from a predefined family of kernels. As an example we consider the family 
$$ \{{g_{\sigma }},\sigma \ge 0\} $$
, where



$$ {g_{\sigma }}(x,y)={\sigma^{-2 }}g\left\{ {\frac{x}{\sigma },\frac{y}{\sigma }} \right\} $$

(8.16)
and the generating kernel 
$$ g $$
is a smoothing kernel [see definition above, (8.13)]. It can easily be seen that for σ > 0, 
$$ {{\hat{g}}_{\sigma }}(0)=1 $$
, that is, 
$$ {g_{\sigma }} $$
is also a smoothing kernel. Besides, 
$$ \mathop{\lim}\limits_{{\sigma \to 0}}F*{g_{\sigma }}(x,y)=F(x,y) $$
, which implies that 
$$ \mathop{\lim}\limits_{{\sigma \to 0}}{g_{\sigma }}=\delta $$
, where 
$$ \delta $$
is a Dirac delta function. These properties allow to construct the image gF defined by



$$ gF(x,y)=\int {F(t,s){g_{{\sigma (x,y)}}}(x-t,y-s)\mathrm{ d}t\mathrm{ d}s}, $$

(8.17)
where we determine 
$$ \sigma $$
depending on current position 
$$ (x,y) $$
. The challenging task is to establish a rule according to which the value 
$$ \sigma $$
can be adapted to the current context of the image. Usually, 
$$ \sigma $$
is required to change as a function of some decisive characteristic such as the gradient, the Laplacian, the curvature, or some other local characteristic of the image. An example of such filtering with the Gaussian function



$$ g(x,y)=\frac{1}{{2\pi }}{{\mathrm{ e}}^{{-\frac{{{x^2}+{y^2}}}{2}}}} $$

(8.18)
as a generating kernel, and the squared deviation 
$$ \xi (x,y) $$
as a local decisive characteristic, is given in Fig. 8.4. Particularly for this example, the current value of 
$$ \sigma $$
was set to



$$ \sigma (x,y)=h\left( {C;\frac{{{\xi_r}(x,y)}}{{{\xi_{\max }}}}} \right), $$

(8.19)
where



$$ h(C;t)=\frac{r}{{1+C{h_1}(t)}} $$

(8.20)
and the function 
$${h_1}:[0,1]\to [0,1]$$
is monotonically increasing from zero to one. For the local squared deviation, we have



$$ \xi_r^2(x,y)=\frac{1}{{2\pi {r^2}}}\int_{{{B_r}(x,y)}} {{{{\left| {F(t,s)-\bar{F}(x,y)} \right|}}^2}\mathrm{ d}t\mathrm{ d}s}, $$

(8.21)




$$ {\xi_{\max }}=\mathop{\max}\limits_{{(x,y)\in \mathbf{D}}}{\xi_r}(x,y), $$

(8.22)
where the ball 
$$ {B_r}(x,y) $$
is a ball of radius 
$$ r $$
located in 
$$ (x,y) $$
; 
$$ \bar{F}(x,y) $$
is the local mean value of 
$$ F $$
inside the ball.

A270527_1_En_8_Fig4_HTML.jpg


Fig. 8.4
Filtering with adaptive bandwidth Gaussian kernel. Left: a patch of the X-ray image shown in Fig. 8.3; middle: the result of the filtering by given tuning parameters 
$$ r $$
and 
$$ C=5r $$
; right: linear filtering with Gaussian kernel by 
$$ \sigma =r $$

The smoothing considered in this example incorporates two tuning parameters, which are the constants 
$$ C $$
and 
$$ r $$
. From (8.19) it follows that 
$$ \sigma $$
varies between 
$$ r $$
and 
$$ r/(C+1) $$
, and the type of this variation depends on the type of the growth of the function h 1 defined in (8.20).

Another class of smoothing techniques is the so-called sigma filtering. Here we use smoothing kernels of the following general form:



$$ {\varphi_{I;x,y }}(t,s)=C(x,y){w_{\tau }}(I(t,s)-I(x,y)){g_{\sigma }}(x-t,y-s), $$

(8.23)
where 
$$ (x,y) $$
is the current update position; 
$$ {g_{\sigma }} $$
and 
$$ {w_{\tau }} $$
are the smoothing kernels, the effective width of which are 
$$ \sigma $$
and 
$$ \tau $$
, respectively; and the local constant 
$$ C(x,y) $$
is determined from the condition



$$ \int\limits_{{{{\mathbb {R}}}^2}} {{\varphi_{{I;x,y}}}(t,s)} \mathrm{ d}t\mathrm{ d}s=1. $$

(8.24)

Because of two kernels acting both in the space and the intensity domains, sigma filtering is sometimes referred to as bilateral filtering (see [4]). As it follows from (8.23), the sigma filter is represented by a family of kernels which depends on two scale parameters 
$$ \sigma $$
and 
$$ \tau $$
. The parameter 
$$ \sigma $$
determines the radius of the ball 
$$ B(x,y) $$
centered in 
$$ (x,y) $$
. Using the parameter 
$$ \tau $$
, it has to be decided which values 
$$ I(x-t,y-s) $$
inside the ball 
$$ B(x,y) $$
have to be weighted down. Normally, these are values which deviate from the value 
$$ I(x,y) $$
. In other words, the kernel 
$$ w $$
generates the local averaging mask that allocates those points within the ball 
$$ B(x,y) $$
where image values are close to 
$$ I(x,y) $$
.

Bilateral filters are especially efficient if applied iteratively. Denoting the filtering operation with 
$$ \rm F $$
, the iterative chain can formally be represented as 
$$ {\rm F_{{{\sigma_n}{\tau_n}}}}{\rm F_{{{\sigma_{n-1 }}{\tau_{n-1 }}}}}\ldots {\rm F_{{{\sigma_1}{\tau_1}}}} $$
.

In [5], it has been shown that the best result by an iterative sigma filter can be achieved if 
$$ {\sigma_1}>\ldots >{\sigma_{n-1 }}>{\sigma_n} $$
” src=”/wp-content/uploads/2016/08/A270527_1_En_8_Chapter_IEq0008104.gif”></SPAN> and <SPAN id=IEq0008105 class=InlineEquation><IMG alt=. It can be shown that such a sequence of filtered images converges to the image which is piecewise constant. As an example, in Fig. 8.5, there are three iterations of bilateral filtering where 
$$ w $$
and 
$$ g $$
are both Gaussian.

A270527_1_En_8_Fig5_HTML.jpg


Fig. 8.5
Sigma filtering. Three iterations of applying sigma filter to the patch of the lung phantom


8.5 Wavelet-Based Nonlinear Filtering


A quite different principle for spatially adaptive noise reduction is based on the reconstruction from selected wavelet coefficients. For the sake of completeness, before describing the related techniques, we will review main points of the theory of the wavelet transform referring mainly to [6] and [7]. In doing so, we will use the terminology of the theory of functions, and therefore, we start by stating formal notations accepted in this theory and used in the text.


$$ {C_m} $$
, 
$$ m\geq 1 $$
, is a space whose elements are 
$$ m $$
times continuously differentiable functions; the space 
$$ {C_0} $$
is the space of continuous functions, and 
$$ {C_{-1 }} $$
, the space of piecewise continuous functions.

Let us have 
$$ j\in \mathbb {Z} $$
and 
$$ m\in \mathbb {N} $$
. As spline space 
$$ S_m^j $$
, we will call the set of all such functions 
$$ f\in {C_{m-2 }} $$
, the restriction of which to the interval 
$$ [{2^j}k,{2^j}(k+1)] $$
is the algebraic polynomial of order 
$$ m-1 $$
; the points 
$$ {2^j}k $$
are called knots of the spline space 
$$ S_m^j $$
.

A space with an inner product 
$$ \left\langle {\cdot, \cdot } \right\rangle $$
is called Hilbert space. In the following, the norm of element 
$$ a\in H $$
is defined by



$$ {{\left\| a \right\|}_H}=\sqrt{{\left\langle {a,a} \right\rangle }}. $$

(8.25)


$$ {L_p}(T) $$
is a space of functions satisfying 
$$ \int_T {{{{\left| {f(x)} \right|}}^p}} <\infty $$
. The space 
$$ {L_2}(T) $$
is a Hilbert space with the inner product



$$ \left\langle {a,b} \right\rangle =\int_T {a(x)b{(x)^{*}}\mathrm{ d}x}, $$

(8.26)
where 
$$ {b^{*}} $$
is a complex conjugate of 
$$ b $$
.


$$ {\ell_p} $$
is a set of all infinite sequences 
$$ {{\{{a_k}\}}_{{k\in \mathbb {Z}}}} $$
such that 
$$ \sum\nolimits_k {{{{\left| {{a_k}} \right|}}^p}<\infty } $$
. The space 
$$ {\ell_2} $$
is a Hilbert space with the inner product



$$ \left\langle {a,b} \right\rangle =\sum\limits_k {{a_k}b_k^{*}}. $$

(8.27)

Span 
$$ {{\{{g_n}\}}_n} $$
is the minimal space spanned on the family 
$$ {{\{{g_n}\}}_n} $$
, i.e., for any 
$$ a\in {\ell_2} $$
,



$$ \sum\nolimits_n {{a_n}{g_n}} =f\in \mathrm{ span}{{\{{g_n}\}}_n}. $$

(8.28)

Operators that map a Hilbert space 
$$ E $$
into Hilbert space 
$$ G $$
will be denoted as 
$$ \mathbf{H}:E\to G $$
. With 
$$ {{\mathbf{H}}^{*}} $$
, we denote the adjoint of 
$$ \mathbf{H} $$
, i.e., such that 
$$ {{\left\langle {\mathbf{H}/f,g} \right\rangle}_G}={{\left\langle {f,{{\mathbf{H}}^{*}}g} \right\rangle}_E} $$
for any 
$$ f\in E $$
and any 
$$ g\in G $$
. The operator is self-adjoint, if 
$$ \mathbf{H}={{\mathbf{H}}^{*}} $$
.

A wavelet is a function 
$$ \psi \in {L_2}(\mathbb {R}) $$
that necessarily satisfies the condition



$$ \int_{\mathbb {R}} {\psi (x)\mathrm{ d}x=\hat{\psi}(0)=0}. $$

(8.29)

Everywhere in the text, it will be supposed that 
$$ \psi $$
is real and normalized, i.e., 
$$ {{\left\| \psi \right\|}_{{{L_2}}}}=1 $$
. Wavelets are “famous” for being well localized both in the space domain and in the frequency domain. The localization properties of wavelets are usually expressed in terms of a so-called space-frequency window. This is a rectangle 
$$ \sigma \times \hat{\sigma} $$
in the space-frequency plane with 
$$ \sigma $$
and 
$$ \hat{\sigma} $$
defined by



$$ {\sigma^2}=\int_{\mathbb {R}} {(x-\bar{x})|\psi (x){|^2}\mathrm{ d}x,} \quad \bar{x}=\int_{\mathbb {R}} {x|\psi (x){|^2}\mathrm{ d}x}, $$

(8.30)




$$ {{\hat{\sigma}}^2}={\pi^{-1 }}\int_0^{\infty } {{{{(\omega -\bar{\omega})}}^2}{{{\left| {\hat{\psi}(\omega )} \right|}}^2}\mathrm{ d}x}, \quad \bar{\omega}={\pi^{-1 }}\int_0^{\infty } {\omega {{{\left| {\hat{\psi}(\omega )} \right|}}^2}\mathrm{ d}x}. $$

(8.31)

The point 
$$ (\bar{x},\bar{\omega}) $$
is the center of the space-frequency window of the wavelet. Since 
$$ \psi $$
is real, the modulus of its Fourier transform is even, and this is why the integration in (8.31) is made over the positive half axis only.

Wavelet 
$$ \psi $$
is called the mother wavelet for the family 
$$ {{\{{\psi_{s,u }}\}}_{{s>0,u\in \mathbb {R}}}} $$
” src=”/wp-content/uploads/2016/08/A270527_1_En_8_Chapter_IEq0008153.gif”></SPAN>, where the function<br />
<DIV id=Equ000832 class=Equation><br />
<DIV class=EquationContent><br />
<DIV class=MediaObject><IMG alt=

(8.32)
is the dilated and translated version of the mother wavelet 
$$ \psi $$
. As an example, let us consider the family generated by so-called Mexican hat wavelet



$$ \lambda (x)=\frac{2}{{\sqrt{3}}}{\pi^{-1/4 }}(1-{x^2}){{\mathrm{ e}}^{{-\frac{{{x^2}}}{2}}}}. $$

(8.33)

In Fig. 8.6 wavelets 
$$ {\lambda_{1,0 }} $$
and 
$$ {\lambda_{5,0 }} $$
of this family (left) as well as their Fourier transforms (right) are shown.

A270527_1_En_8_Fig6_HTML.gif


Fig. 8.6
Mexican hat wavelet. Wavelets 
$$ {\lambda_{s,0 }} $$
for 
$$ s=1,5 $$
(left) and their Fourier transforms (right)

Figure 8.7 shows space-frequency windows of these wavelets. One observes that the bigger is 
$$ s $$
, the wider is the effective width of the wavelet in the space domain and the better is its spectral resolution. In contrary, we obtain better space resolution by smaller 
$$ s $$
. Since 
$$ {{\hat{\sigma}}_s}=\hat{\sigma}/s $$
and 
$$ {\sigma_s}=s\sigma $$
, the area of the space-frequency window does not change.

A270527_1_En_8_Fig7_HTML.jpg


Fig. 8.7
Space-frequency windows of 
$$ {\lambda_{s,0 }} $$
for 
$$ s=1,5 $$

The bivariate function



$$ wf(s,u)=\left\langle {f,{\psi_{s,u }}} \right\rangle,s,u\in \mathbb {R} $$

(8.34)
is called a continuous wavelet transform of the function 
$$ f $$
. Since



$$ \left\langle {\hat{f},{{\hat{\psi}}_{s,u }}} \right\rangle =2\pi \left\langle {f,{\psi_{s,u }}} \right\rangle, $$

(8.35)
the wavelet coefficient 
$$ \left\langle {f,{\psi_{s,u }}} \right\rangle $$
by fixed 
$$ s,u $$
relates to the local contribution made by 
$$ f $$
within the space-frequency window of the wavelet 
$$ {\psi_{s,u }} $$
.

Besides the space-frequency window, another important characteristic of a wavelet is the number of vanishing moments. A wavelet 
$$ \psi $$
is said to have 
$$ m $$
vanishing moments if 
$$ \int {{x^k}\psi (x)\mathrm{ d}x} =0 $$
for 
$$ 0\leq k< m $$
. Since



$$ \frac{{{{\mathrm{ d}}^n}}}{{\mathrm{ d}{\omega^n}}}\hat{\psi}(\omega )=\int {{(-ix)^n}\psi (x){{\mathrm{ e}}^{{-i\omega x}}}\mathrm{ d}x}, $$

(8.36)
the number of vanishing moments of 
$$ \psi $$
is equal to the number of zeros of its Fourier transform at 
$$ \omega =0 $$
. Using this property, it can be shown that a compactly supported wavelet 
$$ \psi $$
has 
$$ m $$
vanishing moments if and only if there exists a compactly supported smoothing function 
$$ \theta $$
such that



$$ \psi (t)={(-1)^n}\frac{{{{\mathrm{ d}}^n}\theta (t)}}{{\mathrm{ d}{t^n}}} $$

(8.37)
(see Chap. 6 of [7] for details); the corresponding wavelet transform of 
$$ f $$
is consequently



$$Wf(s,t) = s^{n+1/2} \frac {{{{\mathrm{ d}}^n}}} {{\mathrm{ d} {t^n}}} (f*{{\bar {\theta}}_s}),$$

(8.38)
where 
$$ {{\bar{\theta}}_s}(t)={\theta_s}(-t) $$
and 
$$\theta_s (t) = s^{-1} \theta (t/s).$$
Hence, the number of vanishing moments is crucial in those applications where the local regularity of the signal has to be measured.

The identity



$$ f(x)=C_{\psi}^{-1}\int_0^{\infty } {\int_{{-\infty}}^{\infty } {\left\langle {f,{\psi_{s,u }}} \right\rangle {\psi_{s,u }}(x)\mathrm{ d}u\frac{{\mathrm{ d}s}}{{{s^2}}}} } $$

(8.39)
with



$$ {C_{\psi }}=\int_0^{\infty } {\frac{{|\hat{\psi}(\omega ){|^2}}}{\omega}\mathrm{ d}\omega } $$

(8.40)
is valid and makes sense for any 
$$ f\in {L_2}\mathbb {R} $$
if 
$$ {C_{\psi }}<\infty $$
. Therefore, the wavelet 
$$ \psi $$
is called admissible if 
$$ {C_{\psi }}<\infty $$
. Condition (8.29) formulated above is necessary for 
$$ \psi $$
to be admissible. However, assuming that 
$$ \psi $$
is “reasonable,” that is, 
$$ \psi $$
is such that it makes sense in the practice, this condition is also sufficient.

There are different possibilities to choose wavelets in 2D. For example, this can be a wavelet 
$$ \psi \in {L_2}({{\mathbb {R}}^2}) $$
that is radially symmetric. Since the wavelet in this case depends effectively on one variable only, the situation is essentially the same as in the 1D case: the 2D wavelet transform is the function that depends on two parameters, the shift (which is now 2D vector) and the dilation. The reconstruction formula is in this case



$$ f{(\bf r)}{\sim}{\mkern6mu} \int_0^{\infty } {\int_{{{{\mathbb {R}}^2}}} {\left\langle {f,{\psi_{s,\bf u }}} \right\rangle {\psi_{s,\bf u }}(\bf r)} \mathrm{ d}\bf u\frac{{\mathrm{ d}s}}{{{s^3}}}}. $$

(8.41)

It is also possible to choose a 2D wavelet that is not radially symmetric. Such wavelets are often called oriented wavelets. The wavelet transform in this case is a function that depends on three parameters, dilation, translation, and rotation, and the corresponding reconstruction formula is



$$ f{(\bf r)}{\sim}{\mkern6mu} \int_0^{{2\pi }} {\int_{{{{\mathbb {R}}}^2}}} {\int_0^{\infty } {\left\langle {f,{\psi_{s,\bf u,\theta }}, } \right\rangle {\psi_{{s,\bf u,\theta }}}(\bf r)\mathrm{ d}\bf u\frac{{\mathrm{ d}s}}{{{s^3}}}\mathrm{ d}\theta } }. $$

(8.42)

In practice, one often applies separable 2D wavelets, i.e., a wavelet of the form 
$$ \varphi (t)\psi (s) $$
, where both 
$$ \varphi $$
and 
$$ \psi $$
are functions from 
$$ {L_2}(\mathbb {R}) $$
, and at least the one of which is a wavelet.

The continuous wavelet transform provides an extremely redundant description of the function 
$$ f $$
. In fact, even a discrete subset of 
$$ {{\{Wf(s,u)\}}_{s,u }} $$
can be used to recover functions exactly. A related discrete wavelet transform is 
$$ W{f_{m,n }}=\left\langle {f,{\psi_{m,n }}} \right\rangle $$
where



$$ {\psi_{m,n }}(x)={a^{-m/2 }}\psi ({a^{-m }}x-nb) $$

(8.43)
for some fixed positive 
$$ a,b $$
.

For a wide range of choices for 
$$ \psi $$
, 
$$ a $$
and 
$$ b $$
discrete families of wavelets constitute frames of the functional space 
$$ {L_2}(\mathbb {R}) $$
. The family 
$$ {{\{{g_n}\}}_n} $$
is called a frame of the Hilbert space 
$$ H $$
if there exist constants 
$$ 0< A\leq B<\infty $$
such that



$$ A\left\| f \right\|_H^2\leq \sum\limits_n {{{{\left| {\left\langle {f,{g_n}} \right\rangle } \right|}}^2}\leq B\left\| f \right\|_H^2} $$

(8.44)
for any 
$$ f\in H $$
. The equivalency 
$$ \left\langle {f,{g_n}} \right\rangle \equiv 0\Leftrightarrow f=0 $$
that follows from this definition implies that frames are complete in 
$$ H $$
, which means that any function in 
$$ H $$
can be approximated by linear combination of frame elements with any degree of accuracy. The frame is called redundant if its elements are linearly dependent.

Only gold members can continue reading. Log In or Register to continue

Stay updated, free articles. Join our Telegram channel

Aug 20, 2016 | Posted by in NUCLEAR MEDICINE | Comments Off on Noise Reduction

Full access? Get Clinical Tree

Get Clinical Tree app for offline access