Image modeling based on a 2-D stochastic subspace system identification algorithm

Abstract

Fitting a causal dynamic model to an image is a fundamental problem in image processing, pattern recognition, and computer vision. In image restoration, for instance, the goal is to recover an estimate of the true image, preferably in the form of a parametric model, given an image that has been degraded by a combination of blur and additive white Gaussian noise. In texture analysis, on the other hand, a model of a particular texture image can serve as a tool for simulating texture patterns. Finally, in image enhancement one computes a model of the true image and the residuals between the image and the modeled image can be interpreted as the result of applying a de-noising filter. There are numerous other applications within the field of image processing that require a causal dynamic model. Such is the case in scene analysis, machined parts inspection, and biometric analysis, to name only a few. There are many types of causal dynamic models that have been proposed in the literature, among which the autoregressive moving average and state-space models (i.e., Kalman filter) are the most commonly used. In this paper we introduce a 2-D stochastic state-space system identification algorithm for fitting a quarter plane causal dynamic Roesser model to an image. The algorithm constructs a causal, recursive, and separable-in-denominator 2-D Kalman filter model. The algorithm is tested with three real images and the quality of the estimated images are assessed.

Appendices

Appendix 1: Kronecker product results

While implementing the 2-D stochastic 4SID algorithm, we will need some Kronecker product identities. In this Section we borrow some results from Magnus and Neudecker (2001) and Martin (2005). Let ${\varvec{A}}\in {\mathbb {R}}^{m \times n}$ and ${\varvec{B}}\in {\mathbb {R}}^{p \times q}$ be two arbitrary matrices. We then have the following results:

$$\begin{aligned} \text {vec}\{{\varvec{A}}^{\top }\}= & {} {\varvec{K}}_{mn} \text {vec}\{{\varvec{A}}\}, \end{aligned}$$

(71)

where ${\varvec{K}}_{mn}$ is called “the commutator matrix” in Magnus and Neudecker (2001) and “the perfect shuffle operator” in Martin (2005), and is a square matrix of dimensions $mn \times mn$, given by

$$\begin{aligned} {\varvec{K}}_{mn}= & {} \left[ \begin{array}{c} I_{mn}(1:m:mn,1:mn) \\ I_{mn}(2:m:mn,1:mn) \\ \vdots \\ I_{mn}(m:m:mn,1:mn) \end{array} \right] . \end{aligned}$$

(72)

${\varvec{K}}_{mn}$ is the matrix obtained from the rows of the $mn \times mn$ identity matrix, $I_{mn}$, by taking every mth row starting with the first, then every mth row starting with the second row, and so on, until the last block obtained by taking every mth row starting with the mth row. Here the MATLAB notation M(a : b : c, d : e) is used to indicate every bth row of M starting with a and ending with c, and columns d through e.

Appendix 2: The operators “vec” and “reshape”

The vec operator converts a matrix into a column vector and the reshape operator, in general, converts a vector into a matrix.

Definition 1

Let ${\varvec{A}}\in {\mathbb {R}}^{m \times n}$. Then $\text {vec}\{{\varvec{A}}\}$ is defined as

$$\begin{aligned} {\varvec{b}}= & {} \text {vec}\{{\varvec{A}}\} \nonumber \\= & {} \left[ \begin{array}{c} {\varvec{A}}(:,1) \\ {\varvec{A}}(:,2) \\ \vdots \\ {\varvec{A}}(:,n) \end{array} \right] \in {\mathbb {R}}^{mn}. \end{aligned}$$

(73)

Definition 2

Let ${\varvec{b}}\in {\mathbb {R}}^{mn}$. Then $\text {reshape}\{{\varvec{b}},m,n\}$ is defined as

$$\begin{aligned} {\varvec{A}}= & {} \text {reshape}\{{\varvec{b}},m,n\} \;=\; \left[ \begin{array}{c|c|c|c} {\varvec{b}}(1:m)&{\varvec{b}}(m+1:2m)&\ldots&{\varvec{b}}((n-1)m+1:mn) \end{array} \right] . \end{aligned}$$

(74)

Note that the definition of $\text {reshape}\{\cdot \}$ in MATLAB does not require ${\varvec{b}}$ to be a vector, but it operates on ${\varvec{b}}$ columnwise.

Appendix 3: Solving the coupled Lyapunov equations

Since (8) is not a Lyapunov equation per se, we need to break it up into its individual components in order to solve it. Adding the symmetry constraints, we obtain a system of equations of the form

$$\begin{aligned} \varPi _h= & {} A_1\varPi _hA_1^{\top } + A_2\varPi _vA_2^{\top } + Q_{hh} \end{aligned}$$

(75a)

$$\begin{aligned} \varPi _v= & {} A_4\varPi _vA_4^{\top } + Q_{vv} \end{aligned}$$

(75b)

$$\begin{aligned} \varPi _{hv}= & {} A_2\varPi _{v}A_4^{\top } + Q_{hv} \end{aligned}$$

(75c)

$$\begin{aligned} \varPi _{h}= & {} \varPi _{h}^{\top } \end{aligned}$$

(75d)

$$\begin{aligned} \varPi _{v}= & {} \varPi _{v}^{\top } \end{aligned}$$

(75e)

$$\begin{aligned} \varPi _{vh}= & {} \varPi ^{\top }_{hv}. \end{aligned}$$

(75f)

Now vectorizing (75a)–(75f), one can find $\text {vec}\{\varPi _{h}\}$, $\text {vec}\{\varPi _{v}\}$, and $\text {vec}\{\varPi _{hv}\}$, by solving the following overdetermined system of equations

$$\begin{aligned} \left[ \begin{array}{c|c|c} I_{n_h^2}-A_1 \otimes A_1 &{} A_2 \otimes A_2 &{} 0_{n_h^2 \times n_hn_v} \\ \hline 0_{n_v^2 \times n_h^2} &{} I_{n_v^2}-A_4 \otimes A_4 &{} 0_{n_v^2 \times n_hn_v} \\ \hline 0_{n_hn_v \times n_h^2} &{} - A_4 \otimes A_2 &{} I_{n_hn_v} \\ \hline I_{n_h^2}-{\varvec{K}}_{n_hn_h} &{} 0_{n_h^2 \times n_v^2} &{} 0_{n_h^2 \times n_hn_v} \\ \hline 0_{n_v^2 \times n_h^2} &{} I_{n_v^2}-{\varvec{K}}_{n_vn_v} &{} 0_{n_v^2 \times n_hn_v} \\ \hline 0_{n_hn_v \times n_h^2} &{} 0_{n_hn_v \times n_v^2} &{} I_{n_hn_v} - {\varvec{K}}_{n_hn_v} \end{array} \right] \left[ \begin{array}{c} \text {vec}\{\varPi _{h}\} \\ \text {vec}\{\varPi _{v}\} \\ \text {vec}\{\varPi _{hv}\} \end{array} \right]= & {} \left[ \begin{array}{c} \text {vec}\{Q_{hh}\} \\ \hline \text {vec}\{Q_{vv}\} \\ \hline \text {vec}\{Q_{hv}\} \\ \hline 0_{n_h^2 \times 1} \\ \hline 0_{n_v^2 \times 1} \\ \hline 0_{n_hn_v \times 1} \end{array} \right] ,\nonumber \\ \end{aligned}$$

(76)

where ${\varvec{K}}_{n_hn_h}\in {\mathbb {R}}^{n_h^2 \times n_h^2}$, ${\varvec{K}}_{n_vn_v}\in {\mathbb {R}}^{n_v^2 \times n_v^2}$, and ${\varvec{K}}_{n_hn_v}\in {\mathbb {R}}^{n_hn_v \times n_hn_v}$ are commutator matrices such that $\text {vec}\{\varPi _{h}^{\top }\}={\varvec{K}}_{n_hn_h}\text {vec}\{\varPi _{h}\}$, $\text {vec}\{\varPi _{v}^{\top }\}={\varvec{K}}_{n_vn_v}\text {vec}\{\varPi _{v}\}$, and $\text {vec}\{\varPi _{hv}^{\top }\}={\varvec{K}}_{n_hn_v}\text {vec}\{\varPi _{hv}\}$.

An alternative solution can be obtained by using a standard discrete Lyapunov equation solver such as “dlyap” from MATLAB. First one would solve

$$\begin{aligned} \varPi _v= & {} \text {dlyap}(A_4,Q_{vv}), \end{aligned}$$

(77)

then solve (75c) for $\varPi _{hv}$, and finally solve

$$\begin{aligned} \varPi _h= & {} \text {dlyap}(A_1,A_2\varPi _vA_2^{\top } + Q_{hh}). \end{aligned}$$

(78)

Once $\{\varPi _h,\varPi _v,\varPi _{hv},S_h,S_v,R\}$ are known, one can compute G and $\varLambda _{0,0}$ from

$$\begin{aligned} G_1= & {} A_1\varPi _h C_1^{\top }+A_2\varPi _v C_2^{\top }+S_h \end{aligned}$$

(79a)

$$\begin{aligned} G_2= & {} A_4\varPi _v C_2^{\top } + S_v \end{aligned}$$

(79b)

$$\begin{aligned} \varLambda _{0,0}= & {} C_1\varPi _h C_1^{\top } + C_2\varPi _v C_2^{\top } + R. \end{aligned}$$

(79c)

Appendix 4: Solving the coupled Riccati equations for $P_h$ and $P_v$

In the 2-D stochastic 4SID algorithm, the state estimate covariance matrices $P_h$ and $P_v$ must be found by solving the coupled Riccati equations (17a)–(17b). To the authors’ knowledge, there is no known algorithm to solve such coupled Riccati equations. An alternative is to solve the two Riccati equations recursively. First, let us compute the state estimation error covariance matrices $\varDelta P_h=\varPi _h-P_h$ and $\varDelta P_v=\varPi _v-P_v$ by direct substitution of the Lyapunov and Riccati equations, i.e.,

$$\begin{aligned} \varDelta P_h= & {} \varPi _h - P_h \;=\; A_1\varPi _hA_1^{\top } + A_2\varPi _vA_2^{\top } + Q_{hh} - A_1P_hA_1^{\top } - A_2P_vA_2^{\top } \nonumber \\&- (G-A_1P_hC_1^{\top }-A_2P_vC_2^{\top })(\varLambda _{00} - C_1P_hC_1^{\top } \nonumber \\&-C_2P_vC_2^{\top })^{-1} (G-A_1P_hC_1^{\top }-A_2P_vC_2^{\top })^{\top } \nonumber \\= & {} A_1\varDelta P_hA_1^{\top } + A_2\varDelta P_vA_2^{\top } + Q_{hh} - (A_1\varPi _hC_1^{\top }+A_2\varPi _hC_2^{\top } + S_h-A_1P_hC_1^{\top }\nonumber \\&-A_2P_vC_2^{\top }) \nonumber \\&\times (C_1\varPi _hC_1^{\top }+C_2\varPi _vC_2^{\top } + R - C_1P_hC_1^{\top } -C_2P_vC_2^{\top })^{-1} (A_1\varPi _hC_1^{\top }+A_2\varPi _hC_2^{\top } \nonumber \\&+ S_h-A_1P_hC_1^{\top }-A_2P_vC_2^{\top })^{\top } \nonumber \\= & {} A_1\varDelta P_hA_1^{\top } + A_2\varDelta P_vA_2^{\top } + Q_{hh} + (A_1\varDelta P_hC_1^{\top } + A_2\varDelta P_vC_2^{\top } + S_h) \nonumber \\&\times (C_1\varDelta P_hC_1^{\top } + C_2\varDelta P_vC_2^{\top } + R)^{-1}(A_1\varDelta P_hC_1^{\top } + A_2\varDelta P_vC_2^{\top } + S_h)^{\top } \nonumber \\ \varDelta P_h= & {} A_1\varDelta P_hA_1^{\top } + \varDelta Q_{hh} + (A_1\varDelta P_hC_1^{\top } + \varDelta S_h)(C_1\varDelta P_hC_1^{\top } + \varDelta R_h)^{-1} (A_1\varDelta P_hC_1^{\top } \nonumber \\&+ \varDelta S_h)^{\top }, \end{aligned}$$

(80)

where

$$\begin{aligned} \varDelta Q_{hh}= & {} A_2\varDelta P_vA_2^{\top } + Q_{hh} \end{aligned}$$

(81a)

$$\begin{aligned} \varDelta S_{h}= & {} A_2\varDelta P_vC_2^{\top } + S_h \end{aligned}$$

(81b)

$$\begin{aligned} \varDelta R_{h}= & {} C_2\varDelta P_vC_2^{\top } + R. \end{aligned}$$

(81c)

Doing the same for $\varDelta P_v$, we obtain

$$\begin{aligned} \varDelta P_v= & {} A_4\varDelta P_vA_4^{\top } + Q_{vv} + (A_4\varDelta P_vC_2^{\top } + S_v) (C_2\varDelta P_vC_2^{\top } + \varDelta R_v)^{-1}(A_4\varDelta P_vC_2^{\top } + S_v)^{\top },\nonumber \\ \end{aligned}$$

(82)

where

$$\begin{aligned} \varDelta R_v= & {} C_1\varDelta P_hC_1^{\top } + R. \end{aligned}$$

(83)

Note that (80)–(82) are now standard algebraic Riccati equations and can be solved with the MATLAB function “dare”. One way to solve (80)–(82) is by first setting $\varDelta P_h=0_{n_h \times n_h}$ and solving (82) to obtain an initial $\varDelta P_v$. Then iterate between (80) and (82) until convergence. The algorithm would be as follows:

Algorithm Riccati

• 1. Initialization: Set $k=0$ and $\varDelta P_h(0)=0_{n_h \times n_h}$.

• 2. Compute $\varDelta P_v(0)$ from (82) with $\varDelta P_h(0)$ as in Step 1.

• 3. Compute

  $$\begin{aligned} \varDelta Q_{hh}(k)= & {} A_2\varDelta P_v(k)A_2^{\top } + Q_{hh} \\ \varDelta S_h(k)= & {} A_2\varDelta P_v(k)C_2^{\top } + S_h \\ \varDelta R_h(k)= & {} C_2\varDelta P_v(k)C_2^{\top } + R \\ \varDelta P_h(k+1)= & {} A_1\varDelta P_h(k)A_1^{\top } + \varDelta Q_{hh}(k) + (A_1\varDelta P_h(k)C_1^{\top } + \varDelta S_h(k)) \\&\times (C_1\varDelta P_h(k)C_1^{\top } + \varDelta R_h(k))^{-1}(A_1\varDelta P_h(k)C_1^{\top } + \varDelta S_h(k))^{\top } \\ \varDelta R_v(k)= & {} C_1\varDelta P_h(k)C_1^{\top } + R \\ \varDelta P_v(k+1)= & {} A_4\varDelta P_v(k)A_4^{\top } + Q_{vv} + (A_4\varDelta P_v(k)C_2^{\top } + S_v) (C_2\varDelta P_v(k)C_2^{\top } \\&+ \varDelta R_v(k))^{-1} (A_4\varDelta P_v(k)C_2^{\top } + S_v)^{\top }. \end{aligned}$$

• 4. Set $k = k + 1$ and return to Step 3 until convergence, i.e., $k \rightarrow k^{\star }$.

• 5. Upon convergence set $\varDelta P_h=\varDelta P_h(k^{\star })$ and $\varDelta P_v=\varDelta P_v(k^{\star })$, where $k^{\star }$ is the convergence iteration index. Compute $P_h$ and $P_v$ from

  $$\begin{aligned} P_h= & {} \varPi _h - \varDelta P_h \end{aligned}$$

  (84a)
  $$\begin{aligned} P_v= & {} \varPi _v - \varDelta P_v. \end{aligned}$$

  (84b)

• 6. Compute the Kalman gain matrices from

  $$\begin{aligned} K_1= & {} (A_1\varDelta P_hC_1^{\top } + A_2\varDelta P_vC_2^{\top } + S_h)(C_1\varDelta P_hC_1^{\top } + C_2\varDelta P_vC_2^{\top } + R)^{-1} \end{aligned}$$

  (85a)
  $$\begin{aligned} K_2= & {} (A_4\varDelta P_vC_2^{\top } + S_v)(C_1\varDelta P_hC_1^{\top } + C_2\varDelta P_vC_2^{\top } + R)^{-1}. \end{aligned}$$

  (85b)

Note that one could use MATLAB’s “dare” function iteratively to find

$$\begin{aligned} \varDelta P_v= & {} \text {dare}(A_4^{\top }, C_2^{\top }, Q_{vv}, \varDelta R_v, S_v) \\ \varDelta P_h= & {} \text {dare}(A_1^{\top }, C_1^{\top }, \varDelta Q_{hh}, \varDelta R_h, \varDelta S_h) \end{aligned}$$

in step 3.