Formal analysis of 2D image processing filters using higher-order logic theorem proving

Two-dimensional (2D) image processing systems are concerned with the processing of the images represented as 2D arrays and are widely used in medicine, transportation and many other autonomous systems. The dynamics of these systems are generally modeled using 2D difference equations, which are mathematically analyzed using the 2D z-transform. It mainly involves a transformation of the difference equations-based models of these systems to their corresponding algebraic equations, mapping the 2D arrays (2D discrete-time signals) over the (z1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_1$$\end{document},z2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_2$$\end{document})-domain. Finally, these (z1\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_1$$\end{document},z2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$z_2$$\end{document})-domain representations are used to analyze various properties of these systems, such as transfer function and stability. Conventional techniques, such as paper-and-pencil proof methods, and computer-based simulation techniques for analyzing these filters cannot assert the accuracy of the analysis due to their inherent limitations like human error proneness, limited computational resources and approximations of the mathematical expressions and results. In this paper, as a complimentary technique, we propose to use formal methods, higher-order logic (HOL) theorem proving, for formally analyzing the image processing filters. These methods can overcome the limitations of the conventional techniques and thus ascertain the accuracy of the analysis. In particular, we formalize the 2D z-transform based on the multivariate theories of calculus using the HOL Light theorem prover. Moreover, we formally analyze a generic (L1,L2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L_1,L_2$$\end{document})-order 2D infinite impulse response image processing filter. We illustrate the practical effectiveness of our proposed approach by formally analyzing a second-order image processing filter.

perform various image processing tasks for controlling the autonomous vehicles, such as noise reduction, color normalization, histogram equalization and edge detection, to enhance the quality of the images captured using various devices such as closedcircuit television (CCTV) and webcams [6]. Similarly, they are widely used in medicine for performing various image pre-and post-processing tasks, such as image quality enhancement, noise removal and image smoothing [5].
The dynamics of these image processing systems are generally modeled using 2D difference equations. Next, the 2D z-transform is used to mathematically analyze these systems. It mainly involves a transformation of the difference equations-based models of these systems to their corresponding algebraic equations, using the definition and various classical properties of the 2D z-transform, while mapping 2D arrays over the ( z 1 ,z 2 )-domain. Finally, these ( z 1 ,z 2 )-domain representations are used to analyze various properties of these image processing systems like transfer function and stability [2].
Conventionally, the image processing filters have been analyzed using paper-and-pencil proof techniques and computer-based symbolic and numerical methods. However, in the former case, the analysis is error-prone due to the highly involved human manipulation, particularly for analyzing the larger and complex image processing systems, and thus we cannot ascertain an absolute accuracy of the analysis in this approach. Similarly, the later approaches suffer from some of their inherent limitations. For example, the symbolic methods involve a large number of unverified symbolic procedures residing in the root of the associated tools [7]. Similarly, the numerical techniques include a finite number of iterations due to the limited power of the computing machines. Moreover, they are based on the mathematical results that are approximated due to the finite precision arithmetic of computers. Therefore, these conventional approaches cannot be trusted when analyzing the image processing filters utilized in various safety-critical areas, such as autonomous driving and medicine.
Formal methods [8] are system analysis techniques that are based on developing a mathematical model of the system using logic and verifying its various properties using deductive reasoning. Higher-order logic (HOL) theorem proving [9,10] is a widely utilized formal method for analyzing many safety-critical systems. In this paper, we propose a HOL theorem proving-based framework for analyzing the image processing filters. In particular, we formalize the 2D z-transform based on the multivariate theories of calculus using the HOL Light theorem prover [11]. The main motivation of selecting HOL Light for the proposed formalization is the presence of the fundamental libraries of multivariate calculus, 1 vectors 2 and matrices, 3 which are required to formally analyze the 2D image processing systems.

Contributions of the paper
The novel contributions of the paper are: • Formalization of 2D z-transform and its region of convergence (ROC). • Formal verification of various classical properties of 2D z-transform, such as linearity, shifting in time-domain, scaling in ( z 1 , z 2 )-domain and complex conjugation. • Formal analysis of a generic ( L 1 , L 2 )-order 2D IIR image processing filter. • Formal analysis of a second-order image processing filter

Preliminaries
This section provides an introduction to the HOL Light theorem prover and the formalization of some fundamental concepts from the multivariate calculus libraries of HOL Light that facilitate the understanding of the rest of the paper.

HOL Light Theorem prover
HOL Light [12] is a proof assistant for developing proofs of the mathematical concepts written as theorems in higher-order logic. HOL Light is implemented in the strongly typed functional programming language ML [13]. A theorem is a statement that is formalized as an axiom or can be implied from already verified theorems using inference rules. Soundness is assured in a theorem proving environment as every new theorem is verified using the primitive inference rules or any other previously verified theorems. HOL Light provides an extensive support of theories, such as Boolean algebra, arithmetic, real numbers, vectors and matrices, which are extensively used in our formalization. Indeed, one of the motivations for selecting the HOL Light theorem prover for the proposed framework is the availability of extensive libraries of vectors and matrices.

Multivariable calculus theories in HOL Light
This section presents an introduction to some fundamental concepts formalized in HOL Light, such as summability, infinite summation and vector summation, and some HOL Light notations that help understanding the rest of the paper.
An N-dimensional vector in HOL Light is formalized as a R N column matrix capturing individual elements as real numbers. All vector operations are then considered as matrix manipulations. Most of the theorems in multivariable calculus theories of HOL Light are proved for functions with an arbitrary data type of R M → R N . Similarly, complex numbers ( C ) can be described as R 2 instead of defining a new data type. The HOL Light symbol &: N → R represents an injection of natural numbers to real numbers. Similarly, the symbol Cx: R → C typecasts real numbers to complex numbers. The symbols Re: C → R and Im: C → R represent the real and imaginary components of a complex number, respectively. The HOL Light symbol % captures the scalar multiplication of a vector or matrix. Similarly, a matrix-vector multiplication is modeled as * * in HOL Light.
The generalized summation over an arbitrary function fn: A → R N is formalized in HOL Light as follows: where vecsum accepts a set st: A → bool over which the summation occurs and a function fn of data type A → R N and returns a generalized vector summation over the set st. Here, the HOL Light function summ provides a finite summation for a fn over real numbers. For example, a mathematical expression We provide the formalization of the summability of a function fn: N → R N over st: N → bool, which ensures that there exist some limit value L: R N , such that ∞ k=0 f (k) = L in HOL Light as: The limit of a function fn: A → R N is formalized as: where the function limt takes a net with components of data type A and a function fn and returns a limit value lt: R N to which fn converges at the given net. It is formalized using the Hilbert choice operator ∈ . Similarly, the concept tends to ( → ) is formalized in HOL Light as: Now, we provide a formalization of an infinite summation, which is used in the formal definition of the 2D z-transform presented in Sect. 5.1.

Definition 6 Infinite Summation of a Function
where the HOL Light function inftsumm accepts st: num → bool specifying the starting point and a function fn of data type N → R N , and returns a limit value lt: R N to which the infinite summation of fn converges from the given st.
Next, we formally verify an equivalence of the infinite summation (Definition 6) to its alternate form in terms of sequential limit as the following HOL Light theorem: Figure 1 depicts our proposed method for analyzing the image processing filters using HOL theorem proving. The user provides the 2D difference equations that model the dynamics of the image processing system, which needs to be analyzed. This 2D difference equation is modeled in higher-order logic using the multivariate calculus theories of HOL Light. In the next step, we formalize the 2D z-transform that is required for mathematically analyzing the image processing systems. It mainly transforms the difference equations-based models of these systems to their corresponding algebraic equations, using the definition and various classical properties, such as linearity, shifting and scaling, of the 2D z-transform, while mapping 2D arrays over the ( z 1 ,z 2 )-domain. Finally, these ( z 1 ,z 2 )-domain representations are used to analyze various properties of these systems, such as transfer function and the solution of the corresponding difference equations.

Formalization of the 2D z-transform
The 2D z-transform of a 2D discrete-time function (2D array) f (n 1 , n 2 ) is mathematically expressed as follows [2]: where f is a function of type N → N → C , and z 1 and z 2 are complex variables. The limits from 0 to ∞ make Eq. (1) as a mathematical representation of a unilateral 2D z-transform. We have opted for this representation based on the same motivation that was considered for one-dimensional z-transform [14] and the Laplace transform [15]. We formalize the 2D z-transform [Eq. (1)] in HOL Light as follows: where z_transform_2d accepts a function of type N → N → C and two complex variables z1: C and z2: C and returns a complex number, which represents the 2D z-transform of f: N → N → C according to Eq. (1).
An essential issue with the applicability of the 2D z-transform of f (n 1 , n 2 ) is the existence of F (z 1 , z 2 ) that occurs due to the presence of the infinite summations in Eq. (1). Thus, we need to identify conditions for the existence of the 2D z-transform. A set of all those values of z 1 and z 2 for which the infinite summations are converging and F (z 1 , z 2 ) is finite (or summable) is known as the ROC. It is mathematically expressed as follows: We formalize the ROC of the 2D z-transform as follows: where ROC_2d accepts a function f: N → N → C and n1 capturing the starting point of the outer summation of the 2D z-transform [Eq. (1)] and returns a set of nonzero values of variables z1 and z2 for which the 2D z-transform of f exists. It is necessary to specify the associated ROC_2d to compute the 2D z-transform. Moreover, the functions z_tr_summable and z_tr_td_summable capture the summability of the function f for the inner and the outer (double) summations, respectively, and are formalized in HOL Light as follows: Moreover, we verify two fundamental properties of ROC, such as the linearity of the ROC and scaling of the ROC, which are quite helpful for formally verifying the classical properties of the 2D z-transform in Sect. 5.2.

Theorem 2 Linearity of ROC
Theorem 2 ensures that if (z1, z2) is inside ROC_2d f n1 and ROC_2d g n1 for functions f and g then it is also inside the intersection of both ROCs for the scaled version of these functions. Similarly, Theorem 3 provides the scaling property with respect to the division by a complex number a.

Formal verification of the classical properties of the 2D z-transform
We use Definitions 7 and 8 and Theorems 2 and 3 for verifying some of the classical properties of the 2D z-transform in HOL Light. This verification plays a vital role in reducing the effort required for analyzing image processing systems, as described later in Sects. 5.3 and 5.4.
Linearity of the 2D z-Transform The linearity of the 2D z-transform is mainly used in decomposing complex (larger) systems to subsystems or combining smaller systems to larger ones having different scaling inputs. It can be mathematically expressed as follows: If Z[f (n1, n2)] = F (z 1 , z 2 ) and Z[g(n 1 , n 2 )] = G(z 1 , z 2 ), then the following holds: The 2D z-transform of a linear combination of 2D sequences (or arrays) is equal to the linear combination of the 2D z-transform of the individual arrays. We verify linearity property in HOL Light as: where a: C and b: C are arbitrary complex constants. Assumptions A1 and A2 capture the regions of the convergence of functions f and g, respectively. The proof of the above theorem is mainly based on Theorem 2 and the linearity of the infinite summation along with some complex arithmetic reasoning. Shifting Property of the 2D z-Transform The shifting property of the 2D z-transform is mostly used for analyzing the 2D linear shift-invariant (LSI) systems. In particular, it is used to solve the difference equations capturing the dynamics of these systems. The shifting property expresses the transform of the shifted signal f (n 1 − m 1 , n 2 − m 2 ) in terms of its 2D z-transform F (z 1 , z 2 ).

Theorem 5 Shifting in Time Domain
where the function in_fst_quad_2d ensures that the function f is nonzero in the first quadrant only and is formalized in a relational form, i.e., f (n1 -m1, n2 -m2), ∀ m1 m2. m1 < n1, m2 < n2. The verification of Theorem 5 is mainly based on the properties of complex numbers along with two properties regarding the negative offset of series and infinite summation. More details about the proof process of this theorem can be found in our proof script. 4 Scaling in (z 1 , z 2 )-domain Property of the 2D z-Transform The scaling property of the 2D z-transform results in shrinking or expansion of the ( z 1 , z 2 )-domain, i.e., 4D complex ( z 1 , z 2 )-plane. If Z[f (n 1 , n 2 )] = F (z 1 , z 2 ) , then two different types of scaling are defined as: If h 1 and h 2 are positive real numbers, then the scaling is interpreted as expansion of the 4D complex ( z 1 , z 2 )-plane. On the other hand, multiplication by w 1 −n 1 and w 2 −n 2 [Eq. (6)] shrinks the ( z 1 , z 2 )-domain. We verify the above theorems in HOL Light as: Complex Conjugation Property of the 2D z-Transform The complex conjugation property facilitates an easy manipulation of the 2D z-transform of conjugated arrays. It is mathematically expressed as follows: where f * (n 1 , n 2 ) represents the complex conjugate of an array f (n 1 , n 2 ) . The corresponding formalization of the complex conjugation property in HOL Light is given as follows:

Formal verification of a ( L 1 , L 2 )-order 2D infinite impulse response (IIR) image
processing filter 2D digital filters [1] are integral components of the image processing systems. Their main responsibility includes the decomposition of an image to multiple frequency bands, restricting a 2D array/signal to a certain frequency band and providing the input-output relationship of these systems. For example, a low-pass filter allows a range of frequencies less than a certain threshold [2]. The analysis of an image processing filter mainly involves developing its mathematical model using a 2D difference equation. The next step is to apply 2D z-transform on both sides of the difference equation. Finally, the definition and the classical properties of the 2D z-transform are used to perform transfer function-based analysis of the given filter.
The impulse response of a discrete-time system captures its behavior for the scenario when dirac-delta function is acting as an input array [2]. 2D image processing infinite impulse response (IIR) filters have a nonzero impulse response function over an infinite length of time. For these filters, the present output depends on the present input and all previously computed input and output values.
Mathematically, the 2D image processing filters are described using the following difference equation [16]: where α(l 1 , l 2 ) and β(k 1 , k 2 ) are input and output coefficients, respectively. The output array y(n 1 , n 2 ) is a linear combination of the previous K 1 − 1 and K 2 − 1 output samples, the present input x(n 1 , n 2 )) , and L 1 − 1 and L 2 − 1 previous input samples. Moreover, for the shift-invariant filter, α(l 1 , l 2 ) and β(k 1 , k 2 ) are the complex constants ( C ). Therefore, Eq. (8) is known as a linear constant coefficient difference equation (LCCDE). The 2D z-transform of a (L 1 , L 2 ) th difference represented in the form of f (n 1 , n 2 ) is given as: The corresponding transfer function of the 2D IIR filter is mathematically expressed as [16]: To formally verify the transfer function of the 2D filter [Eq. (10)], we formalize the (L 1 , L 2 ) th difference as follows: The function l1l2th_difference accepts a function f: N → N → C , coefficients of the difference equation c l1 l2, the order (L1, L2) of the 2D difference equation and the variables n1 and n2 and returns the (L 1 , L 2 ) th difference. It uses the function vsum s f twice to capture the double summation.
Next, we formalize a general LCCDE [Eq. (8)] as follows: Next, we verify the 2D z-transform of the (L 1 , L 2 ) th difference [Eq. (9)] as: The 2D z-Transform of the (L 1 , L 2 ) th Difference (8) y(n 1 , n 2 ) = where Assumption A1 ensures that (z 1 , z 2 ) are in the region of convergence of the function f. Assumption A2 implies that the function f is in the first quadrant. Finally, the conclusion provides the 2D z-transform of the (L 1 , L 2 ) th difference. The verification of the above theorem is mainly based on induction on N1 and N2 and Theorems 2 and 4 along with the following lemma about the summability of (L 1 , L 2 ) th difference equation.
To verify the transfer function of the 2D filter [Eq. (10)], we have to ensure that the 2D input and output arrays exist in the first quadrant only. Moreover, the denominator of Eq. (10) should be nonzero. We formalize both these requirements in HOL Light as follows: Assumption A1 provides the ROC for LCCDE. Assumption A2 ensures that the input and output 2D arrays are in the first quadrant. Assumption A3 captures the time-domain model of the 2D IIR filter, i.e., the LCCDE (Eq. (8)). Finally, the conclusion presents the transfer function of the 2D IIR filter. The proof process of the above theorem is based on the linearity and shifting properties of the 2D z-transform (Theorems 4 and 5) and summability of the (L 1 , L 2 ) th difference (Lemma 1) along with some complex arithmetic reasoning. Theorem 10 provides the transfer function of a generic 2D IIR image processing filter and is quite useful in the verification of the second-order 2D medical image processing filter described in Sect. 5.4.

Formal verification of a second-order 2D image processing filter
To illustrate the practical utilization and effectiveness of the proposed formalization of the 2D z-transform, we apply it to formally analyze a second-order image processing filter that is widely used for performing various tasks, such as noise removal [1], image smoothing [2] and quality enhancement [5].
A second-order image processing filter is graphically represented by the flow graph shown in Fig. 2. A flow graph is a collection of branches (directed connections) and nodes (input and output 2D arrays), where nodes are connected using branches. The constants c 01 , c 10 , c 11 , c 02 , c 12 , c 20 , c 21 and c 22 in Fig. 2 represent the gains of each branches, whereas z 1 −1 and z 2 −1 present the shift right (horizontal delay) and shift up (vertical delay) operations, respectively. We can mathematically describe this filter using the following linear difference equation.