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.

Similar content being viewed by others
Explore related subjects
Discover the latest articles, news and stories from top researchers in related subjects.Notes
It is well known that images are nonstationary processes, however, under mild violation of these assumptions, the 2-D 4SID algorithm performed very well with real images, as seen in Sect. 7.
In order to avoid an abuse of notation when there is no time index involved, and for lack of better words, we still use the terms past and future as is commonly used in time domain subspace system identification.
For zero mean processes, \({\varvec{A}}|{\varvec{B}}={\varvec{A}}{\varvec{D}}{\varvec{B}}^{\top }({\varvec{B}}{\varvec{D}}{\varvec{B}}^{\top })^{-1}{\varvec{B}}\) denotes the orthogonal projection of \({\varvec{A}}\) onto \({\varvec{B}}\), where \({\varvec{D}}\) is a diagonal weight matrix.
4SID stands for State Space Subspace System IDentification.
There is no need to compute the Q part, just the triangular factor L will suffice, thus achieving a significant reduction in the number of computations.
M(1 : m, : ) denotes the first m rows of M and all its columns.
MATLAB is a registered trademark of The Math Works, Inc.
References
Abramatic, J., Attasi, S., Chieze, J., & Curien, N. (1975). Mod‘eles statistiques bidimensionnels et application au traitement num’erique des images. In Colloque national sur le traitement du signal et ses applications, Nice, France.
Aggarwal, G., Chowdhury, A., & Chellappa, A. (2004). A system identification approach for video-based face recognition. In: Proceedings of the international conference on pattern recognition, Washington, D.C.
Alata, O., & Olivier, C. (2003). Choice of a 2-D causal autoregressive texture model using information criteria. Pattern Recognition Letters, 24, 1191–1201.
Ali, M., Chughtai, S. S., & Werner, H. (2009). Identification of spatially interconnected systems. In Proceedings of the 48th IEEE conference on decision and control, held jointly with the 2009 28th Chinese control conference (pp. 7163–7168).
Attasi, S. (1976). Modelling and recursive estimation for double indexed sequences. In R. Mehra & D. Lainiotis (Eds.), System identification: Advances and case studies (pp. 289–348). New York: Academic Press.
Azimi-Sadjadi, M., & Wong, P. (1987). Two-dimensional block Kalman filtering for image restoration. IEEE Transactions on Acoustics, Speech, and Signal Processing, 35, 1736–1749.
Barry, P., Gran, R., & Waters, C. (1976). Two-dimensional filtering: A state space approach. In Proceedings of the IEEE conference on decision and control, Clearwater, FL.
Bose, N. (1977). Problems and progress in multidimensional system theory. Proceedings of the IEEE, 65, 824–840.
Camcioglu, E., & Kayran, A. H. (2013). Spectrum estimation using new efficient 2-d AR lattice modeling of random fields. In Africon, 2013 (pp. 1–6).
Chen, C., Tsai, J., & Shieh, L. (2003). Modeling of variable coefficient Roesser’s model for systems described by second order partial differential equation. Circuits, Systems, and Signal Processing, 22, 423–463.
Delp, E., Kashyap, L., & Mitchell, O. (1979). Image data compression using autoregressive time series models. Pattern Recognition, 11, 313–323.
Dillabough, M. (2007). Discrete state-space models for distributed parameter systems (Unpublished master’s thesis). Sudbury, ON: School of Graduate Studies Laurentian University.
Ding, T., Sznaier, M., & Camps, O. (2006). Robust identification of 2-D periodic systems with applications to texture synthesis and classification. In Proceedings of the IEEE conference on decision and control, San Diego, CA.
Fornasini, E., & Marchesini, G. (1978). Doubly-indexed dynamical systems: State-space models and structural properties. Mathematical Systems Theory, 12, 59–72.
Fraanje, R., & Verhaegen, M. (2005). A spatial canonical approach to multidimensional state-space identification for distributed parameter systems. In Proceedings of the international workshop on multidi-mensional systems, Wuppertal.
Galkowski, K. (2000). Higher order discretization of 2-d systems. IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications, 47(5), 713–722. doi:10.1109/81.847876.
Givone, D., & Roesser, R. (1972). Multidimensional linear iterative circuits—General properties. IEEE Transactions on Computers, 21, 1067–1073.
Golub, G., & Van Loan, C. (1996). Matrix computations (3rd ed.). Baltimore, MD: John Hopkins University Press.
Habibi, A. (1972). Two-dimensional Bayesian estimate of images. Proceedings of the IEEE, 60, 878–883.
Hunger, R. (2007). Floating point operations in matrix-vector calculus (Technical Report No. Version 1.3). Technische Universität München Associate Institute for Signal Processing. Retrieved from https://mediatum.ub.tum.de/doc/625604/625604.pdf.
Jain, A. (1981). Advances in mathematical models for image processing. Proceedings of the IEEE, 69, 502–528.
Kaczorek, T. (1985). Two dimensional linear systems. Berlin: Springer.
Katayama, T. (2005). Subspace methods for system identification. Berlin: Springer.
Katayama, T., & Kosaka, M. (1979). Recursive filtering algorithm for a two-dimensional system. IEEE Transactions of Automatic Control, 24, 130–132.
Katayama, T., & Picci, G. (1999). Realization of stochastic systems with exogenous inputs and subspace identification methods. Automatica, 35, 1635–1652.
Kayran, A. H., & Camcioglu, E. (2014). New efficient 2-d lattice structures for general autoregressive modeling of random fields. IEEE Transactions on Signal Processing, 62(6), 1590–1602.
Kizilkaya, A., & Kayran, A. (2005). Estimation of 2-D ARMA model parameters by using equivalent AR approach. Journal of the Franklin Institute, 342, 39–67.
Kuduvalli, G., & Rangayyan, R. (1993). An algorithm for direct computation of 2-D linear prediction coefficients. IEEE Transactions on Signal Processing, 41, 996–1000.
Kung, S. Y., Levy, B. C., Morf, M., & Kailath, T. (1977). New results in 2-D system theory, part II: 2-D state-space models-realization and the notions of controllability, observability, and minimality. Proceedings of the IEEE, 65, 945–961.
Kurek, (1985). The general state-space model for a two-dimensional linear digital system. IEEE Transactions on Automatic Control, 30(6), 600–602.
Kurek, & Zaremba, M. (1993). Iterative learning control synthesis based on 2-D system theory. IEEE Transactions on Automatic Control, 38, 121–125.
Lashgari, B. (1981). Two-dimensional approximation, model reduction and minimum variance error estimation. Unpublished doctoral dissertation. University of Southern California, Los Angeles, CA.
Lashgari, B., & Silverman, L. (1986). Recursive estimation of two-dimensional processes with application to image processing. Circuits Systems and Signal Processing, 5, 212–226.
Lashgari, B., Silverman, L., & Abramatic, J. (1983). Approximation of 2-D separable in denominator filters. IEEE Transactions on Circuits and Systems, 30, 107–121.
Lele, S., & Mendel, J. (1987). Modeling and recursive state estimation for two-dimensional noncausal filters with applications in image restoration. IEEE Transactions on Circuits and Systems, 34, 1507–1517.
Lin, T., Kawamata, M., & Higuchi, T. (1987). Decomposition of 2-D separable-denominator systems: Existence, uniqueness, and applications. IEEE Transactions on Circuits and Systems, 34, 292–296.
Magnus, J., & Neudecker, H. (2001). Matrix differential calculus with applications in statistics and econometrics. New York: Wiley.
Mardia, K. V., Kent, J. T., & Bibby, J. M. (1979). Multivariate analysis. New York: Academic Press.
Martin, C. (2005). Higher-order Kronecker products and tensor decompositions. Unpublished doctoral dissertation, Cornell University, Ithaca, NY.
Moonen, M., De Moor, B., Vandenberghe, L., & Vandewalle, J. (1989). On and off line identification of linear state space models. International Journal of Control, 49, 219–232.
Owens, D., Amann, N., Rogers, E., & French, M. (2000). Analysis of linear iterative learning control schemes-a 2-D systems/repetitive processes approach. Multidimensional Systems and Signal Processing, 11, 125–177.
Ramos, J. (1994). A subspace algorithm for identifying 2-D separable-in-denominator filters. IEEE Transactions on Circuits and Systems, 41, 63–67.
Ramos, J., Alenany, A., Shang, H., & Lopes dos Santos, P. (2011). Subspace algorithms for identifying separable in denominator 2-D systems with deterministic inputs. IET Control Theory and Applications, 5, 1748–1768.
Ramos, J., & dos Santos, P. (2011). Subspace system identification of separable-in-denominator 2-D stochastic systems. In Proceedings of the IEEE conference on decision and control and European control conference, Orlando, FL.
Rao, C. R. (1964). The use and interpretation of principal component analysis in applied research. Sankhya, A, 26, 329359.
Roesser, R. (1975). A discrete-state-space model for linear image processing. IEEE Transactions on Automatic Control, 20, 1–10.
Strintzis, M. (1976). Comments on “two-dimensional Bayesian estimate of images”. Proceedings of the IEEE, 64, 1255–1257.
Suresh, B., & Shenoi, B. (1981). New results in two-dimensional Kalman filtering with applications to image restoration. IEEE Transactions on Circuits and Systems, 28, 307–319.
USC-SIPI Image Database, (n.d.). http://sipi.usc.edu/database.
Van Overschee, P., & De Moor, B. (1996). Subspace identification for linear systems: Theory, implementation, and applications. New York: Kluwer.
Verhaegen, M. (1994). Identification of the deterministic part of mimo state space model given in innovations form from input–output data. Automatica, 30(1), 61–74.
Verhaegen, M., & Verdult, V. (2007). Filtering and system identification: A least squares approach. Cambridge: Cambridge University Press.
Wang, Z., Bovik, A., Sheikh, H., & Simoncelli, E. P. (2004). Image quality assessment: From error visibility to structural similarity. IEEE Transactions on Image Processing, 13, 600–612.
Woods, J., & Ingle, V. (1981). Kalman filtering in two dimensions: Further results. IEEE Transactions on Acoustic, Speech, and Signal Processing, 29, 188–197.
Wu, Z. (1985). Multidimensional state-space model Kalman filtering with application to image restoration. IEEE Transactions on Acoustic, Speech, and Signal Processing, 33, 1576–1592.
Xiao, C., Sreeram, V., Liu, W., & Venetsanopolos, A. (1998). Identification and model reduction of 2-D systems via the extended impulse response gramians. Automatica, 34, 93–101.
Yang, R., Xie, L., & Zhang, C. (2006). H2 and mixed H2/H1 control of two-dimensional systems in Roesser model. Automatica, 42, 1507–1514.
Zhang, J. Y., & Steenaart, W. (1990). A very fast kalman filter for image restoration. In Proceedings of the IEEE international symposium on circuits and systems, New Orleans, LA.
Author information
Authors and Affiliations
Corresponding author
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:
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
\({\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
Definition 2
Let \({\varvec{b}}\in {\mathbb {R}}^{mn}\). Then \(\text {reshape}\{{\varvec{b}},m,n\}\) is defined as
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
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
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
then solve (75c) for \(\varPi _{hv}\), and finally solve
Once \(\{\varPi _h,\varPi _v,\varPi _{hv},S_h,S_v,R\}\) are known, one can compute G and \(\varLambda _{0,0}\) from
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.,
where
Doing the same for \(\varDelta P_v\), we obtain
where
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
in step 3.
Rights and permissions
About this article
Cite this article
Ramos, J.A., Mercère, G. Image modeling based on a 2-D stochastic subspace system identification algorithm. Multidim Syst Sign Process 28, 1133–1165 (2017). https://doi.org/10.1007/s11045-016-0385-4
Received:
Revised:
Accepted:
Published:
Issue Date:
DOI: https://doi.org/10.1007/s11045-016-0385-4