Entropic Image Restoration as a Dynamic System with Entropy Operator

Entropy and the classical variational principle of the statistical physics are the effective tools for modeling and solving a lot of applied problems. There are many definitions of "entropy" functions. The book by Kapur (1989) contains some of them. The classical definition of the physical entropy was introduced by L.Boltzmann Boltzmann (1871) and was developed for Fermiand Einstein-statistics Landau & Livshitz (1964). Notion of entropy was introduced for para-statistics that have position between Fermiand Einstein-statistics (Ohnuki & Kamefuchi (1982), Dorofeev et al. (2008)).


Introduction
Entropy and the classical variational principle of the statistical physics are the effective tools for modeling and solving a lot of applied problems. There are many definitions of "entropy" functions. The book by Kapur (1989) contains some of them. The classical definition of the physical entropy was introduced by L. Boltzmann Boltzmann (1871) and was developed for Fermi-and Einstein-statistics Landau & Livshitz (1964). Notion of entropy was introduced for para-statistics that have position between Fermi-and Einstein-statistics (Ohnuki & Kamefuchi (1982), Dorofeev et al. (2008)).
Variation principle of entropy maximization turned out very useful for information theory, the base of which connected with Shannon (Shannon (1948), Kullback & Leibler (1951)). This direction is developed in the book , where introduced generalized information entropies by Fermi-Dirac and Bose-Einstein (entropy with parameters). Entropy maximization are applied to image reconstruction from projections Byrne (1993). A large number of applications of the entropy maximization principle is contained in Fang et al. (1997), Maslov (2003).
In these papers the entropy conditional maximization problems with linear constraints equalities were considered only. However there are many problems of entropy maximization with feasible set that is described by a system of inequalities and not only a linear one.
In this paper we design the models of the entropy image reconstruction from projections (EIRP) as the entropy linear (ELP) and quadratic maximization problems (EQP), where the feasible sets are described by the system of the equalities and inequalities of appropriate types (linear and quadratic one).
The regular procedure for design of multiplicative algorithms with p-active variables with respect to dual variables and to mixed (dual and primal) variables proposed for the problem solving. The choice of the active variables is implemented by feedback control with respect to the current state of the iterative process.
The problem of reconstruction of images of the objects distorted by noises and hidden from direct observation arises in the different fields. One of the trends in the solution of the problem is based on the tomographic investigation of an object, i.e., the construction of its layer-by-layer projections. The projections can be formed as external irradiation sources(X-ray, ultrasonic sources) and internal ones (positron emission) as also with the aid of their combination (nuclear magnetic resonance) (Herman (1980), Dhawan (2003)). In Popkov (1997) is shown that a distribution of the absobed photons in slab maximizes the generalized entropy by Fermi-Dirac under the set of projections. A generalization consists in the inclusion of an additional parameters in entropy function, through which it is possible to take into account prior information on the object.
Our contribution to the theory and applications of the EIRP consists in three parts.
The first contribution is the general entropy models in the terms of entropy linear programming (ELP) or entropy quadratic programming (EQP) that underlie in the static procedures of computer tomography. At the beginning of the static procedure, it is occured an accumulation of a complete set of the projections by means of the external irradiation of the object. Then it is solved the ELP or EQP. As a result, for the prescribed prior image, we obtain an entropy-optimal restored image, which we will be called a posterior image.
It calls for a rather high irradiation intensity so as to afford a sufficient noise immunity of a reconstructed image. However, for some classes of tomographic investigations a high irradiation intensity is extremely undesirable.
Multiplicative procedures represent to the ELP and EQP solving. Apparently, the first general approach for synthesis of such procedures was proposed in (Dubov et al. (1983)). Simple multiplicative algorithm was applied to minimization of strictly convex functions on nonnegative orthant. Later, the multiplicative algorithms with respect to dual variables are used for solving conditional minimization and mathematical programming problems (Aliev et al. (1985), Popkov (1988), Popkov (1995a). Also, the multiplicative algorithms are used for solving nonlinear equations (Popkov (1996)). The multiplicative procedures for finding nonnegative solutions of the minimization problems over nonnegative optant were proposed again in the paper (Iusem et al. (1996)).
Some types of the multiplicative algorithms are derived from approach based on the Bregman function and generalized projections with Shannon's entropy. In this case we obtain so-called row-action algorithms, iterations of which have a multiplicative form. The algorithms of this type was developed and modified (Herman (1982), Censor (1981), Censor (1987), Byrne (1996), Censor & Zenios (1997)).
It is necessary to note that in the most cited works the multiplicative algorithms are applied to the problems of entropy maximization with linear constraints equalities. We consider the ELP problem, where a feasible set is described by the system of the linear equalities and inequalities. The regular procedure for design of multiplicative algorithms with p-active variables is proposed for solving of this problem. On the basis of the procedure above we sinthesize the algorithms with respect to dual variables and to mixed (dual and primal) variables simultaneously. The choice of the active variables is implemented by feedback control with respect to the current state of the iterative process. Convergence study of the multiplicative algorithms is based on the continuous analogues of the algorithms and equivalence of the iterative sequences generated by the dual and mixed type algorithms. (Popkov (2006)).
Our second contribution is connected with a basically another approach to the EIRP. It is a dynamic procedure consisting in the sequential refinement in time of the image synthesized. The suggested procedures do not require a high irradiation intensity and display a high noise stability.
Entropic Image Restoration as a Dynamic System with Entropy Operator 3 On the each step t of the dynamic procedure t-posterior image is build up as a solution of the ELP or EQP, using the current t-prior image and the current projection. The (t − s), (t − s + 1),...,t-posterior images take part in formation of (t+1)-prior image. So the dynamic procedures are procedures with feedback.
The dynamic procedures are closed in the sense that at each stage for the current t-prior image and the t-projection, the entropy-optimal t-posterior image is built up, by which the (t+1)-prior image is corrected.
We consider a diverse structures of the dynamic procedures with feedback and investigate their properties. The example of application of these procedures is presented.
It is shown that the proposed dynamic procedures of the EIRP represent the dynamic systems with entropy operator (DSEO). And our third contribution is an elements of the qualitative analysis of the DSEO. We consider the properties of the entropy operator (boundedness, Lipschitz constant).

Mathematical model of the static EIRP procedure
Consider a common diagram of monochrome tomographic investigation ( fig. 1), where external beams of photons S irradiate the flat object in the direction AB. The object is monochromatic, and is described by the two-dimensional function of optical density ψ(x, y) in the system of Cartesian coordinates. Positive values of the density function are limited: (2.1) The intensity of irradiation (projection) w at the point B of the detector D ( fig. 1) is related by the Radon transformation: where the integration is realized along straight AB.
The tomographic procedure form some feasible sets for the vectorψ: where L(ψ) is the h-vector function, and g is the h-vector. We consider the quadratic approximation of the function L: where L is the (h × m)-matrix with nonnegative elements l ki ≥ 0; Q(ψ) is the h-vector of the quadratic forms: where Q k is the symmetric (m × m)-matrix with elements q k ij ≥ 0.

47
Entropic Image Restoration as a Dynamic System with Entropy Operator Now it is returned to the projection function (2.2), and we use its quadratic approximation: where: w(B)={w(B 1 ),...,w(B n )}, T is the (n × m)-matrix with elements t ki ≥ 0; F(ψ) is the n-vector function with components F rψ )=ψ ′ F rψ , where F r is the symmetric (m × m)-matrix with elements f r ij ≥ 0. Any tomographic investigation occurs in the presence of noises. So the n-projections vector u is a random vector with independent components u n , n = 1, m, where u 0 is the ideal projections vector (without noise), and σ 2 is the dispersion of the noise. It is assumed that the dispersions of the noise components are equal.
Thus, the feasible set D(ψ) is described the following expressions: -the projections are Tψ + F(ψ)=u, (2.7) -the possible set of the density vectors is The class P of the density vectors is characterized by the following inequalities: It is assumed that among dimensions of the density vectors (m), the projection vectors (n), and the possible set (h) the following inequality exists: It is assumed that the feasible set is nonempty for the class (2.9), and there exists a set of the density vectorsψ (2.9), that belong to the feasible set D (2.7, 2.8).
We will use the variation principle of the EIRP Popkov (1997), according to which the realizable density vector (function)ψ maximizes the entropy (the generalized information entropy by Fermi-Dirac): where: -E = {E 1 ,...,E m } is the m-vector characterizing the prior image (prior probabilities of photon absorption in the object); If there is information about more or less "grey" object then we can use the next entropy function (the generalized information entropy by Boltzmann): where e = 2, 73.
Thus, the problem of the EIRP can be formulated in the next form: where the feasible set D(ψ) is described by the expressions (2.7 -2.8) and the class P is described by the inequalities (2.9). This problem is related to the EQP or the ELP depend on the feasible set construction.

Statements and algorithms for the ELP and the EQP
Transform the problem (2.13) to the general form, for that introduce the following designations: Then the problem (2.13) takes a form: under the following constraints: Remark that the constraints (2.9) are absent in the problem (3.2, 3.4), as they are included to the goal function.

The ELP problem
1. Optimality conditions. The feasible set in the ELP problem is described by the next expressions: (3.7) Consider the Lagrange function for the ELP (3.2, 3.6): whereλ,μ are the Lagrange multipliers for constraints-equalities and -inequalities (3.6) correspondingly. Assume that the Slater conditions are valid, i.e., there exists a vector x 0 such that Lx 0 <ĝ, Tx 0 =û.
According to Polyak (1987) the following expressions give the necessary and sufficient conditions optimality of the triple (x,λ,μ) for the problem (3.2 -3.5): where ⊗ designates a coordinate-wise multiplication.
The following designations are used in these expressions: From the optimality conditions (3.9) we have: The Lagrange multipliersμ and the exponential Lagrange multipliers z = exp(λ) are defined by the next equations and inequalities: 2. Multiplicative algorithms with (p+q)-active variables An active variables are vary at the sth iteration, and the remaining variables are not vary. We will consider multiplicative algorithms with respect to dual variables (z,μ) for solution of the system (3.16). At the each step of iteration it will be used p components of the vector z, and q components of the vector μ. The number of the active variables is valid to the next relation: The multiplicative algorithms with p+q-active variables can be represented in the following form: The parameters γ, α are the step coefficients. In Popkov (2006) the multiplicative algorithms in respect to the mixed type (prime and dual variables) are introduced, and the method of the convergence of these algorithms are proposed.

The EQL problem
1. Optimality condition. Consider the EQL problem (3.2 -3.5), and introduce the Lagrange function: According to the optimality conditions (3.9, 3.10) we have: Transform these equations and inequalities to the conventional form in which all variables are nonnegative one: x l q k jl .
(3.23) 2. Multiplicative algorithms of the mixed type with (p+q+w)-active variables. We use p active prime x variables, q active dual variables z for the constraints-equalities, and w active dualμ variables for the the constraints-inequalities. The algorithm takes a form: The parameters β, γ, α are the step coefficients.

Active variables.
To choice active variables we use feedback control with respect to the residuals on the each step of iteration. Consider the choosing rule of the active variables for the ELP problem (3.2, 3.3, 3.4). Introduce the residuals One of the possible rules is a choice with respect the maximum residual. In this case it is necessary to select p maximum residual ϑ i 1 ,...,ϑ i p and q maximum residual ε k 1 ,...,ε k q for the each iterative step s. The numbers i 1 ,...,i p and k 1 ,...,k q belong to the intervals [1, n] and [1, h] respectively.
Now we define the rule of the (p+q)-maximal residual in the following form: According to this rule we have the chain of inequalities: We can see that all dual variables are sequentially transformed to active ones during I + J + 2 iterations It is repeated with a period of I + J + 2.

Dynamic EIRP procedure with feedback
The basic idea of the dynamic procedure lies in the sequential (stage by stage) obtaining of the projections and the solution of the sequence of the appropriate ELP or EQP. The feasible sets in these problems consist on two subsets. One of them describes the class of the possible density functions (2.7), and it is not depends from irradiation of the object. The other subset depends on the measured projections (2.8).
In the general case the procedure holds where L is the feedback operator, which characterizes the transformation of the t-prior image, and t, (t − 1),..., (t − s)-posterior images to the (t + 1)-prior image.
Represent the operatorL in the following form: where ǫ is a small positive real number.
Then the dynamic EIRP procedure (the discrete procedure) takes a form: * ) ,...,ψ ((t−s), * ) ). (4.4) Now let the variable t is continuous one. Then under ǫ → 0 we will have the continuous dynamic EIRP that is described be the differential equation: dE(t) dt = L(E t ,ψ (t, * ) ,...,ψ ((t−s), * ) ). (4.5) The t, (t − 1),..., (t − s)-posterior images are defined by the ELP or EQP problems, which represent the appropriate entropy operators (4.1). So the dynamic EIRP procedure (4.1, 4.2, (4.4)) represents the discrete dynamic system with entropy operator (the discrete DSEO) and its the continuous analog (4.5) represents the continuous dynamic system with entropy operator (the continuous DSEO). Some general properties of the DSEO will be described in the next section.

Structures of the dynamic EIRP procedures
Let us consider some partial cases. One of them relates to the examination of a Markov version of the procedure when the information only on the t-posterior image is used to shape up E (t+1) : E (t+1) = L(E t ,ψ t, * ). In the second case, information collections at t, t − 1, . . . , t − s stages are used for the estimation of the current meanψ t of the posterior image: Finally, in the third case, information collections at t, t − 1, . . . , t − s stages are used for the estimation of the current meanψ t and dispersion d t of the posterior image: (4.8) We will introduce the following types of the dynamic procedures of the EIRP: • the identical feedback (I − f eedback) (4.9) • the feedback with respect to the current mean of image (CM − f eedback) (4.11) • the feedback with respect to the current mean and dispersion of image (CMD − f eedback) (4.14)

Investigation of the dynamic EIRP procedure with I-feedback
Consider the problem (2.13) in which the feasible set is the polyhedron, a = 0, b = 1, and the constraints to the possible density functions (2.7) are absent. In this case t-posterior density function hold: The exponential Lagrange multipliers z 1 ,...,z n are defined from the following equations: According to the definition of the I-feedback procedure we have: The iterative process (4.17) can be considered as the method of the simple iteration applying to the eguations: Then the system of equations (4.19) has the unique solution E * .
If Ψ ′ (0)=1/a ≥ 1, (a ≤ 1), then the auxiliary equation has the unique solution, and the method of the simple iteration is converged to this solution.
Now it is necessary to find a conditions when ϕ i [z(E)] ≤ 1. The sufficient conditions for it is formed by the following theorem.
Theorem 2. Let the matrix T in (2.7) has the complete rank n, and the following conditions be valid: Proof. Consider the Jacobian of the vector-functionΦ(z t ). Its elements take a form: The equality to zero is reached when z → ∞. So, the functions Φ 1 ,...,Φ n are strictly monotone decreasing.

57
Entropic Image Restoration as a Dynamic System with Entropy Operator www.intechopen.com

Computer experiment
Consider the dynamic EIRP when the tomographic device gives the orthogonal linear projections with the matrix As a test image, use is made of the LENA test (IEEE Image Processing), on which the spot is placed (fig. 2, the left upper window).
To the right and below window, the projections with noise are shown. In the example, the noise/signal ratio amounted to 0.3.
It is necessary to restored the LENA with the spot having the noisy projections. We use "the pure LENA", which is shown in the second upper window, as the 0-prior image E 0 = {E 0 1 ,...,E 0 m }. This problem of the EIRP is described by the ELP with the constrains-equalities. The multiplicative algorithms with 1-active dual variable (3.18) is used for solution of the problem. The modification of the dynamic procedure with I-feedback involved the following. At each stage t we shaped up the auxiliary vector In parallel, the current mean of the components of the vectorψ t, * are calculated: The modified dynamic procedure takes the form: In fig. 2, in the right upper window, the result of the EIRP with the static procedure is shown. In the middle lower window, fig. 2 shows results of the EIRP by the dynamic procedure with I -feedback, and in the right lower window results of the modified dynamic procedure is shown. The quality of the right image is obvious.
The test image 2 is shown in the fig. 3.

Dynamic systems with entropy operator (DSEO)
We can see that dynamic procedures of the image restoration from projections represent a dynamic discrete system with the particular type of the entropy operator with the generalized entropy Fermi-Dirac (3.2), and the feasible set that is described by the inequality and the equality (3.4).

59
Entropic Image Restoration as a Dynamic System with Entropy Operator

www.intechopen.com
In general case the class of the continuous DSEO is described by the following differential equations:

y[t, u(t), v(t)])
, with the entropy operator: where: H(t, y | u) is an entropy function with the H-parameters u; the vectors (y, u) ∈ R n , v ∈ R m , and D(t, v) is a feasible set depended from the D-parameters v.
In these equations U is the n-vector-function, and V is the m-vector-function.

Classification of the DSEO
Some physical analogues we will use for construction of the classificatory graph. In particular, from the equations (5.1, 5.2) it is seen the rates of parameters is proportional to the flows. The entropy function is a probability characteristics of a stochastic process. So the H-parameters are the parameters of this process.
We will use the following classificatory indicators: The classificatory graph is shown in the fig. 4. At the beginning we consider some properties of the entropy operator, notably, the HD, B, Eq, Plh -entropy operator that is included to the HD -DSEO: where Boltzmann-entropy function is and The (n × m)-matrixT =[t ki ≥ 0] has a full rank equal n.  The HD, B, Eq, Plh -entropy operator describes the mapping of the sets U m + (u − , u + ) and V n + (v − , v + ) into the set Y⊂R m + of the operator's values. We will characterize this mapping by two local Lipschitz-constants -L U and L V , i.e. (5.9) We will use the upper estimations of local Lipschitz-constant: where Y U and Y V are the U-Jacobian and the V-Jacobian of the operator y[u, v] respectively.
Evaluate the normalized entropy operator in the following form: The matrix T in (5.11) has a full rank n and the normalized elements, i.e. ∑ n k=1 t ki = 1 for all i = 1, m. Also it is assumed that the condition of the dominating diagonal is valid for the quadratic matrix TT ′ , i.e. the following inequalities take a form: (5.14)

61
Entropic Image Restoration as a Dynamic System with Entropy Operator www.intechopen.com The feasible set D = {x : T x = v, x ≥ 0} is not empty, notable, there exists some subset of the nonnegative vectors x ∈D.
Designate the Lipschitz-constant for the normalized operator (5.11) asL U andL V respectively, i.e.
We will use the upper estimations of local Lipschitz-constant: (5.16) where X U and X V are the U-Jacobian and the V-Jacobian of the operator x[u, v] respectively.
Thus we have the following equalities: Thus, we will calculate the local Lipschitz-constants estimations for the normalized entropy operator (5.11, 5.12) and then apply the formulas (5.17, 5.18).
The normalized entropy operator x(u, v) can be represented by the form: m, (5.19) where the Lagrange multipliers λ j (u, v), (j = 1, n) as the implicit functions from u, v define by the equations: (5.20) 1. Estimations of the norm's matrix X U . The (m × m)-matrix X U takes a form: We will use Euclidean vector norm ( y 2 ), with which two matrix norm are consisted (see Voevodin (1984)): -the spectral norm where σ max is the maximal eigenvalue of the matrix A; It is known that It is assumed that X U = X U 2 . We have from (5.19) the following equality: (5.23) and the n × m-matrix m (5.24) In these expressions ⊗ is coordinate-wise multiplication of the vector's components to the rows of the matrix.
According to (5.21) and the relation between the spectral norm and the Euclidean norm, we have: (5.26) x max = max (i,u,v) x Now consider the equations (5.20), and differentiate the left and right sides of these equations by u. We obtain the following matrix equation: From this implies that Here the (n × n)-matrix Φ λ has elements According to (5.30) we have Thus the norm's estimation of the matrix X U takes a form:

Estimations of the norm's matrix
It is assumed that X V = X V 2 . We have from (5.19) the following equality: where the (m × n)-matrix X λ = −x ⊗ T ′ ; (5.36) and the n × n-matrix According to (5.36) and the relation between the spectral norm and the Euclidean norm, we have: Now consider the equations (5.20), and differentiate the left and right sides of these equations by v. We obtain the following matrix equation: Here the (n × n)-matrix Φ λ is defined by (5.31). According to (5.41) we have Thus the norm's estimation of the matrix X V takes a form: So, we can see that it is necessary to construct the norm's estimation for the matrix Φ −1 λ ,a s well as for the norm's estimation of the matrix X U .
3. Estimation of the spectral norm of the matrix Φ −1 λ . The matrix Φ λ (5.31) is symmetric and strictly negative defined for all λ. Therefore, it has n real, various, and negative eigenvalues (see Wilkinson (1970)). We will order them in the following way: The spectral norm of an inverse matrix is equal to the inverse value of the modulus of the maximum eigenvalue μ max of the initial matrix, i.e.
To definite the value M we resort to the Gershgorin theorem (see Wilkinson (1970)). According to the theorem any eigenvalue of a symmetric strictly negative definite matrix lies at least in one of the intervals with center −c k (λ) and the width 2 ρ k (λ): We can apply the lower estimation for the left side of this interval using (5.14): where x min = min (i,u,v) x i (u, v). (5.50) Thus, in a view of (5.28), we have (5.51)

Estimation of the local Lipschitz-constants.
According to (5.34) and (5.51), the estimation of the local U-Lipshitz-constant for the normalized entropy operator (5.19) takes a form: where the vector z max has the components z max .
For realization of the this way it is necessary to construct the majorant of the operator (5.61, 5.62). We use the following inequality Bellman (1961): Then for the operator (5.62) the following estimate is valid: where the matrix C has the elements t ki t ji , (k, j)=1, n.
(5.68) It is follows from (5.68) that the matrix C takes a form: Thus, we can consider in the capacity ofẑ the nonnegative solution of the equation: The general solution of the equation (5.70) can be represented in the following form:

Conclusions
Many applied problems can be formulated as the ELP or EQP, models of which it is proposed in the paper. The multiplicative algorithms with p-active variables and feedback are the effective methods of their solution. The dynamic procedure of the image restoration from projections (IRP) increase appreciably the quality of the restored image in the presence of noise in the measurements. It is represented a classification of the dynamic procedures and it is investigated a stability of the procedure with I-feedback. Also it is shown that in general case the dynamic procedure of the IRP is the dynamic system with entropy operator (EO). The analytical-numerical methods investigation the problem of the EO-boundedness and calculation of the Lipschitz constant for EO are proposed.