2D Iterative MAP Detection: Principles and Applications in Image Restoration

618 D. KEKRT, T. LUKE ˇ S, M. KL´ IMA, K. FLIEGEL, 2D ITERATIVE MAP DETECTION 2D Iterative MAP Detection: Principles and Applications in Image Restoration Daniel KEKRT, Tom´ aˇ s LUKE ˇ S, Miloˇ s KL´ IMA, Karel FLIEGEL Dept. of Radioelectronics, Czech Technical University in Prague, Technick´ a 2, 166 27 Prague 6, Czech Republic {kekrtd1, lukestom, klima, fliegek}@fel.cvut.cz Abstract. The paper provides a theoretical framework for the two-dimensional iterative maximum a posteriori detection. This generalization is based on the concept of detection algorithms BCJR and SOVA, i.e., the classical (one-dimensional) iterative detectors used in telecommunication applications. We generalize the one-dimensional detection problem considering the spatial ISI kernel as a twodimensional finite state machine (2D FSM) representing a network of the spatially concatenated elements. The cellular structure topology defines the design of the 2D Iterative decoding network, where each cell is a general combinationmarginalization statistical element (SISO module) exchanging discrete probability density functions (information metrics) with neighboring cells. In this paper, we statistically analyse the performance of various topologies with respect to their application in the field of image restoration. The iterative detection algorithm was applied on the task of binarization of images taken from a CCD camera. The reconstruction includes suppression of the defocus caused by the lens, CCD sensor noise suppression and interpolation (demosaicing). The simulations prove that the algorithm provides satisfactory results even in the case of an input image that is under-sampled due to the Bayer mask. Keywords Iterative detection, 2D iterative decoding netwoks, maximum a posteriori probability criterion, defocus suppression, deconvolution, denoising, de-mosaicing, binary image restoration, image processing. 1. Introduction The iterative detection based on the criterion of maximum a posteriori probability (MAP) is broadly used within the Turbo code detection [10]-[13]. The principle of Turbo decoder is closely related to the concept of the Viterbi algorithm [18]. Specifically, it includes the forward-backward recurrent algorithm (FBA-BCJR [16]), or the Viterbi algorithm with soft output (SOVA [17]) for code detection with binary input data. These systems work with one-dimensional signals where the independent variable represents time. The basic idea of two-dimensional iterative data detection was presented in [9]. Thiennviboon et al. applied iterative detection to digital image half-toning [7], [8]. The efficiency of the iterative detection has been shown on simulations with black and white defocused images [1], [3] and motion blurred images [2], where Kekrt et al. also introduced new variants and concatenations of network topologies. The aim of the paper is to provide a theoretical framework for the iterative detection in two-dimensional form with emphasis on applications in the field of image processing. We extended the layered topology and provide a comparison of the existing topologies with respect to their applicability for image reconstruction. The input data are considered to be two-dimensional and they are coded by a spatial kernel. The spatial kernel introduces a general inter-symbol interference (ISI) to the signal. The spatial ISI channel can be seen as a two-dimensional finite state machine (2D FSM) and it can be decomposed into a horizontal-vertical network of simple combinational logic (cells). There are several possible ways how a coding network can be arranged. The topology of the coding network is further a template for an iterative decoding network where each cell is created by a general combination-marginalization statistical element or a SISO module. The goal of such a detector is to reconstruct the original (input) data, thus the inter-symbol interference and additive noise are suppressed. The detection process is given by a sequential activation of all cells of the network in a predefined order. Within this process, the cells are exchanging soft information (the discrete probability density) between each other. In the first part of the paper, the theoretical framework of the two-dimensional iterative detection is described. Further various topologies of detection network are described including their activation schedules. An application for the reconstruction of black and white images is shown using selected topologies. The algorithm performs deconvolution, noise suppression and interpolation (de-mosaicing) of the binary image. The performance of the network topologies is statistically analysed by Monte Carlo method using random data with different levels of noise and blur. At the end, computational complexity and implementation issues are discussed. RADIOENGINEERING, VOL. 23, NO. 2, JUNE 2014 619 2. General System Model We assume a two-dimensional system (Fig. 1). The matrix Dcontains discrete values of input data. In the ISI channel, generally two-dimensional, inter-symbol interference is introduced into the input data. The output of the ISI channel is denoted by the matrix Q. The elements of the matrix Qare in discrete values. They are further referred to as code symbols. The Qmatrix is modified in a front-end channel block defined by the function fs(.). This block provides an unambiguous mapping S=fs(Q)from code symbols to the signal that is transmitted to the transmission channel. The front-end channel block transformation function fs(.)is memory-less and often linear. In the transmission channel, the signal Sis corrupted by a signal interference W. This stochastic process can be described as an additive uncorrelated noise. In the model, we assume additive white Gaussian noise (AWGN). The corrupted signal Renters into the detector, which is composed from a network front-end (a gateway) and an iterative decoding network (IDN). The IDN works with probability densities instead of isolated values. Therefore, there is a gateway in front of the IDN that calculates these probability densities (a posteriori metrics) based on the obtained realization Rand a known noise distribution in the transmission channel. The metrics are further denoted by S(ˇ ξ)because it may refer to probability or its version transformed by an unambiguous mapping. Index Fmarks metrics considered as a posteriori (“forward”) metrics. Index Bmarks a priori (“backward”) metrics. An argument of each soft metric consists of a testing estimator ˇ ξ. It represents a possible value from a range of values of a given random variable and it is used to address a given metric. 2.1 Deterministic ISI Channel The two-dimensional ISI channel can be defined by the convolution: q[k,l] = f(G,N[k,l]) (1) =∑ k0,l0∈L(G) g[k0,l0]d[k+k0,l+l0] where G={g[k0,l0]}k0,l0∈L(G)is a convolution kernel of the channel, L(G)denotes the cardinality of the kernel and N[k,l] = {d[k+k0,l+l0]}k0,l0∈L(G)is the convolution region of input data d[k,l]∈Adwhere kand lare spatial indexes in both dimensions of the image space. The input data are assumed to be binary (Ad={0,1}). The output of the channel q[k,l]represents data interference within the region N[k,l] and it depends on the kernel size and type. Under the assumption of discrete input data, q[k,l]∈Aqis also discrete in values. All possible values that may result from the convolution create an alphabet Aq(G). 2.2 Cellular Models of ISI Channels The convolution (1) is discrete operation due to the discrete inputs and it can be realized by a two-dimensional net- / Noises W data Source Q SB(ˇ Q) D ˆ D R SF(ˇ Q) Estimated S Detector data channel (2D FSM) Channel front-end fs(.) ISI (gateway) front-end Network decoding network Iterative channel Random Fig. 1. Block diagram of a general two-dimensional transmission system. / d[k, l] l−1 R[k, l] C[k, l] k k−1 k−2 ll −2 / d[k, l] B[k, l] R[k, l] C[k, l] l−2 k k−1 k−2 ll −1 (a) Fixed topology. (b) Extended fixed topology. Fig. 2. The shapes of state variables R[k,l],C[k,l]and B[k,l] for a general reduced ISI channel G(R) 3×3. work of mutually connected elements. The network forms the finite state machine (2D FSM). These elements are shown in Fig. 3a. Let us assume that they have a general set of inputs and outputs: VI|O[`]∈{V(i) I|O[`]}i∈L(VI|O[`]) with a discrete range of values V(i) I|O[`]and cardinality L(VI|O[`]). Values of all inputs and outputs of one element create a set NIO =[ ` VI|O[`](2) and input values of one element are denoted as N=[ ` VI[`].(3) The output values of the element {VO[`]}`=f(G,N)(4) are a function of the kernel Gand the inputs N. Firstly, the elements of the network can be understood as associated combinatorial logic that forms a Mealy type machine. Therefore, Fig. 3a contains more output ports although the output of the ISI channel (1) is only one. These “extra” outputs can be seen as state variables. They also represent input states of the neighboring elements. Secondly, the elements create a combinatorial logic where the state variables are not present. Then the element will have one output and the cardinality card(NIO) = card(N)+1. There are three basic network topologies that are able to perform the convolution (1). These topologies are called: fixed topology, variable topology and layered topology. The fixed topology (FT) is based on a connection of combination elements in horizontal and vertical direction [11]. An extended version (EFT) has also diagonal connections [1]. In both cases, the topology does 620 D. KEKRT, T. LUKE ˇ S, M. KL´ IMA, K. FLIEGEL, 2D ITERATIVE MAP DETECTION not depend on the shape of the kernel Gand NIO[k,l] = {R[k,l],C[k,l],B[k,l],d[k,l],R[k,l+1],C[k+1,l],B[k+ 1,l+1],q[k,l]}. The first half of the set is given by the inputs and the second half are the outputs. The elements are exchanging state variables R,C, and Bin vertical, horizontal and diagonal direction. The state variables contain a part of the convolution region N[k,l]and may be mutually disjoint (but it is not required). The initial states R[k,l], C[k,l], and B[k,l]together with the data d[k,l]create the whole convolution region N[k,l]. The input data are situated in the right down corner of the region. The final states R[k,l+1],C[k+1,l],B[k+1,l+1]⊂N[k,l], transmitted to the neighboring elements, are given by a shift of information in the indicated direction (Fig. 2a). The state variables have to be placed and shifted in such a way that the network meets the condition of causality. Fig. 2b shows an example of the state variables for the reduced 3×3 kernel (i.e. a kernel with zero values in the corners). The variable topology (VT) [11] represents a different approach to the solution of (1). The VT network is dependent on the shape of the kernel. It contains two types of combination elements (so-called broadcasters and receivers) and generally higher number of interconnections. The processing is carried out in two levels. On the lower level, the input data are copied into the state variables c[k,l+l0,`] = d[k,l] by broadcasters. The broadcasters (with the input-output set N(B) IO [k,l] = d[k,l]∪{c[k,l,`]}) create all mutually overlapped convolution regions N[k,l] = {c[k+k0,l+l0,`]}. Consequently, these regions are combined into the outputs using elements on the higher level (receivers). The set of input-output variables of the receivers is given as NIO[k,l] = N[k,l]∪q[k,l]. The FT and VT topologies are able to model the convolution kernel of an arbitrary shape. On the contrary, the layered topology can be applied only if the convolution (1) is separable into two orthogonal directions. The following condition must hold true for this type of kernels G=G(V)×G(H), where ×denotes the Cartesian product. Separable kernel always has a square or rectangular shape and the convolution (1) can be decomposed into a pair of one-dimensional convolutions c[k,l] = fH(G(H),N(H)[k,l]) = ∑l0∈L(G(H))gH[l0]d[k,l+l0]and q[k,l] = fV(G(V),N(V)[k,l]) = ∑k0∈L(G(V))gV[k0]c[k+k0,l], where N(H)[k,l] = {d[k,l+l0]}l0∈L(G(H))and N(V)[k,l] = {c[k+k0,l]}k0∈L(G(V))are one-dimensional convolution regions. Although the applicability of the layered topology is limited, there are several advantages, which we will discuss later in the paper. Both layers of the model can be implemented using two layers with variable topology or two layers with fixed topology. In the latter case, inputs and outputs of the lower layer create a set N(H) IO [k,l] = {R[k,l],d[k,l],R[k,l+1],c[k,l]} and for the upper layer N(V) IO [k,l] = {C[k,l],c[k,l],C[k+ 1,l],q[k,l]}. 2.3 Random IECS-ML Channel For the signal Stransmitted through a memory-less (ML) channel with independent eliminated channel states (IECSML), we may write R=S+W(5) where Wrepresents an additive noise. Values of noise are assumed to be spatially invariant in the space [k,l]. The inputoutput relation of the random channel may be rewritten in the following scalar form r[k,l] = s[k,l]+w[k,l](6) where each noise value w[k,l]has the same probability distribution pw(ξ). The probability distribution is considered to be the Gaussian with non-zero mean value. 3. 2D MAP Detection The two-dimensional maximum a posteriori (MAP) detection is in principle similar as in the one-dimensional case. More about the one-dimensional MAP detection can be found in [11, 12, 13]. We will focus on a transition of the problem from the single-stage (optimal) MAP detection to the iterative (suboptimal) MAP detection, which can be realized in practice. 3.1 Single-Stage (Optimal) MAP Detection The optimal 2D MAP detector is based on the following MAP criterion ˆ d[k,l] = arg M  ˇ d[k,l] M  ˇ D:ˇ d[k,l] S(R,ˇ D)!(7) where ˆ d[k,l]is the output estimate of the reconstructed signal at the [k,l]position, ˇ d[k,l]denotes a testing estimator (a possible value of reconstructed signal), ˇ D:ˇ d[k,l]is a set of all possible realizations, which contains at the given position ˇ d[k,l]and M as well as M represents general marginalization operators. Marginalization process runs through joint soft measures of the detector S(R,ˇ D)and quantifies the degree of truth of the statement ˇ D=Dwhen a realization Rwas received. The amount of metrics S(R,ˇ D)grows exponentially with the size of ˇ D. The basis of this growth is the cardinality of the input alphabet. Therefore, the MAP criterion can not be implemented directly. In order to make the problem feasible, it is necessary to factorize the criterion into smaller parts. It is possible due to the assumption that the noise values w[k,l]in the matrix Ware not correlated. It stands that pW(Ξ) = ∏k,lpw(ξ)and the metric of the detector S(R,ˇ D) = C  k,lSF(r[k,l]|ˇ N[k,l])C C  k,lSB(ˇ d[k,l]) =SF(R|ˇ D)C SB(ˇ D)(8) can be decomposed into individual, overlapping, convolution regions ˇ N[k,l], where SF(r[k,l]|ˇ N[k,l]) are a posteriori RADIOENGINEERING, VOL. 23, NO. 2, JUNE 2014 621 metrics and SB(ˇ d[k,l]) are a priori metrics of data. We assume that the data are spatially invariant [k,l]. If the data has uniform distribution, the component SB(ˇ D)in (8) can be neglected and the MAP criterion becomes the Maximum Likelihood (ML) criterion. All factorized metrics in (8) are combined into the joint metric of the detector using a general combination operator C . Combination and marginalization operators can be divided according to an operation domain and a detection technique. It is the symbol detection (SyD) or the sequential detection (PgD). The symbol detection tries to minimize the error detection of a symbol d[k,l]. The sequential (or also page) detection tries to minimize the error detection of all data Das a whole. Tab. 1 shows an overview of the operators where M(.) = −lnP(.) and min∗(x,y) = min(x,y)−ln(1+e−|x−y|). Md-PgD version is numerically the most effective because it has the lowest numerical complexity of the operators. Domain Detect. M M C C C −1 Probability Page max max Π× ÷ Probability Symbol max Σ Π × ÷ Metric Page min min Σ+− Metric Symbol min min∗Σ+− Tab. 1. Summary of combination and marginalization operators. 3.2 Iterative (Suboptimal) MAP Detection It has been shown that the MAP criterion can be for ML-IECS channel decomposed into individual convolution regions ˇ N[k,l]. This leads to the iterative detection, which basically reduces the complexity of the single-stage MAP criterion from the size of the input data Dto the size of the kernel G. The iterative detector (decoding network) is a system of mutually connected elements within one of the topologies described in Sec. 2.2. Each element of the network decodes the particular convolution region. The minimum number of network elements is equal to the size of the input data D. Combination-marginalization elements are shown in Fig. 3b. Hereafter, we will refer to them as SISO (Soft-In Soft-Out) modules. In essence, they represent statistical soft inversion in relation to the elements of the ISI channel model (Fig. 3a). Each port of the SISO module is bidirectional and carries the whole discrete soft metric (probability density) quantifying the reliability of all possible values of the variable (not only one value of the variable as in the case of the ISI channel model). During the detection process, each module initially take over the input metrics SIfrom its neighbors and combine them into partial joint a posteriori metrics S(ˇ NIO). The joint a posteriori metrics S(ˇ NIO)denote reliabilities of all possible arrangements (realizations) of the particular convolution region. Subsequently, the SISO module marginalize the set of joint metrics into outputs. This operation is called an activation of the SISO module. SISO modules are activated progressively according to the so-called activation schedule. The activation schedule has to have a logic that corresponds to the information flow (state and auxiliary variables) within the ISI channel model. If the activation schedule is selected inappropriately (against the flow of the state variables), the time of the convergence to the correct solution (given in the number of iterations) increases. The term iteration will be understood as a time interval in which each SISO module is activated at least once. Each iteration is concluded with the exchange of metrics SO(ˇ VI|O[`]) →SI(ˇ VI|O[`]) (9) SI(ˇ VI|O[`]) ←SO(ˇ VI|O[`]) between SISO modules. / VI[1] VI[0] VO[L] VO[L−1] VO[`+ 1] VI(`) f(.) (a) General processing element. / SI(ˇ VI[1]) SO(ˇ VI[1]) SI(ˇ VI[0]) SO(ˇ VI[0]) SO(ˇ VO[L]) SI(ˇ VO[L]) SI(ˇ VO[L−1]) SO(ˇ VO[L−1]) SI(ˇ VO[`+ 1]) SO(ˇ VO[`+ 1])SI(ˇ VI[`]) SO(ˇ VI[`]) f−1(.) (b) Soft-In Soft-Out module. Fig. 3. A general processing element and its soft inversion (SISO module). At the beginning of the detection process, the input metrics SIare set according to the a priori knowledge about the variables which they represent. If the a priori knowledge is unknown, these metrics are set uniformly (each possible value of the testing estimator has the same probability). In this case, only forward (input) metrics SFfrom the front-end of the iterative detector will be non-uniform. The iterative detection is suboptimal. The final estimate can be regarded as optimal after infinite number of iterations. In most cases, however, the convergence of the network is fast (in order of units of iterations). When the amount of information exchanged between SISO modules stabilizes, it does not change much with the following iterations. The system reaches a state when the metrics are “locked” in the feedback (iterative loops). The detection process is then terminated and the final hard decision is made ˆ VI|O[`] = arg M  ˇ VI|O[`]SI(ˇ VI|O[`]) C SO(ˇ VI|O[`])(10) by a combination of the input and output metrics on the particular port of the SISO module (10). 622 D. KEKRT, T. LUKE ˇ S, M. KL´ IMA, K. FLIEGEL, 2D ITERATIVE MAP DETECTION 4. Iterative Decoding Networks In Sec. 3.2, we have described the base of the iterative detection, cellular structure of the two-dimensional iterative MAP detector and its functional blocks. The activation of the SISO module will be further discussed in detail together with the integration of the SISO module within different topologies (described in Sec. 2.2). 4.1 Network Elements - SISO Modules The activation of the SISO module is composed of two steps - combination and marginalization. The combination merges input (a priori and a posteriori) soft metrics {SI(ˇ VI|O[`])}`of all input-output variables VI|O[`]into a joint a posteriori metric S(ˇ NIO)using a particular combination operator C . This operation is given as S(ˇ NI|O) = C  ˇ VI|O[`]∈ˇ NIO SI(ˇ VI|O[`]) (11) where the condition ˇ VI|O[`]∈ˇ NIO arises from the combination table of the SISO module. A posteriori metrics S(ˇ NIO)∈ {S(N(i) IO )}icorrespond to the likelihood of individual realizations of the set of inputs and outputs N(i) IO . The alphabet of all realizations of the set NIO is denoted as AIO ={N(i) IO }i. For example, a 3 ×3 channel with a binary input has card(AIO) = 512. Therefore, a set of 512 joint metrics is obtained in the first step of the activation of the SISO module. i7→AIO ic[0]7→Acic[2]ic[1]iq7→Aq 0 0 0 0 0 1 1 0 0 1 2 0 1 0 3 3 1 1 0 4 4 0 0 1 1 5 1 0 1 2 6 0 1 1 4 7 1 1 1 5 Tab. 2. Table of combinations controlling the example simple SISO module. Joint metrics are marginalized into the output (a posteriori) metrics {SO(ˇ VI|O[`])}`by the marginalization operator M . This second step of the activation is defined as SO(ˇ VI|O[`]) =  M  ˇ NIO:ˇ VI|O[`] S(ˇ NIO) C −1SI(ˇ VI|O[`]) (12) where ˇ NIO :ˇ VI|O[`]⊂AIO denotes a set of all inputs and outputs that contain the testing estimator ˇ VI|O[`]. Subsequently, the SISO module sends the output metrics into the neighboring modules. The mechanism is explained in the following example. We assume a simple one-dimensional kernel G(H)= {g0g g0}. The output of such a channel is q=g0(c[0]+ c[1])+gc[2], where c[0],c[1], and c[2]are binary inputs. For testing estimators of inputs, outputs, and input-output sets, we may write ˇc[`]∈Ac={0,1}, ˇq∈Aq={0,g0,2g0,g,g+ g0,g+2g0}and ˇ NIO ={ˇc[0],ˇc[2],ˇc[1],ˇq}. Individual variables are associated into combinations (Tab. 2) where i are indexes of mappings into alphabet of the given variable. Each out of eight table rows represents one realization N(i) IO of the set ˇ NIO. Four densities are the input of the combination. Three two-element SI(ˇc[`]) and one sixelement SI(ˇq). The combination of input metrics gives eight joint metrics S(ˇ NIO). The joint metrics are marginalized into three two-element metrics SO(ˇc[`]) and one six-element metric SO(ˇq). An example of the first step of the activation process is calculation of the joint metric S(ˇ NIO = {1,1,0,g+g0}) = SI(ˇc[0] = 1)C SI(ˇc[2] = 1)C SI(ˇc[1] = 0)C SI(ˇq=g+g0). This combination is marked by blue squares in the table (Tab. 2). An example of the second step of the activation process can be the calculation of the output metric SO(ˇc[0] = 0)=(S(ˇ NIO ={0,0,0,0})M S(ˇ NIO = {0,1,0,g})M S(ˇ NIO ={0,0,1,g0})M S(ˇ NIO ={0,1,1,g+ g0})) C −1SI(ˇc[0] = 0). This marginalization is marked by red squares. Each other SISO module uses the same principle. A difference is only in input-output variables, their cardinality and combination table that controls the combinationmarginalization process. 4.2 Fixed Topologies Fig. 4a. shows the iterative decoding network with fixed topology. The topology does not depend on the shape of the kernel and is composed of the SISO modules depicted in Fig. 4b. Basic module connections are horizontal and vertical [11]. An extended version has also diagonal connections [1]. The activation schedule of the detector is in Fig. 4a (the purple dashed line). The activation can be performed row by row or column by column or in “zig-zag” manner. All these variants progress in the same direction as the flow of the state variables in the ISI channel model. Networks with fixed topology are also called as networks with the marginalization at the symbol block level. The convolution region is separated into blocks. The advantage of the fixed topology is high efficiency due to state variables that often have a significant overlap. Large amout of information has to be transmitted and stored between the modules. Optimal separation of the convolution region into the state variables is difficult for some kernel types. Sometimes it is not even possible. In that case, the state variables have a higher cardinality and the memory requirements are raising. The extended variant with diagonal connections tries to suppress these problems. The state variable then allows higher variability for the convolution region decomposition. There are three state variables, so they can have lower cardinality and the memory requirements are reduced in exchange for slight degradation of the performance. RADIOENGINEERING, VOL. 23, NO. 2, JUNE 2014 623 / l−1l l + 1 k+ 1 k k−1 f−1(.)f−1(.) f−1(.)f−1(.)f−1(.) f−1(.)f−1(.) f−1(.) f−1(.) schedule Activation (a) The IDN topology with marked activation schedule. / SF(ˇ d)[k, l] SB(ˇ d)[k, l] SI(ˇ R)[k, l] SO(ˇ R)[k, l] SF(ˇq)[k, l] SB(ˇq)[k, l] SI(ˇ R)[k, l + 1] SO(ˇ R)[k, l + 1] SO(ˇ C)[k, l] SO(ˇ C)[k+ 1, l] SI(ˇ C)[k+ 1, l] SI(ˇ B)[k+ 1, l + 1] SO(ˇ B)[k+ 1, l + 1] SI(ˇ B)[k, l] SO(ˇ B)[k, l] SI(ˇ C)[k, l] f−1(.) (b) The IDN cell in the node [k,l]. Fig. 4. Iterative decoding network marginalizing at the symbol block level. 4.3 Variable Topologies The variable topology can not be depicted in a general form because it depends on the shape of the kernel. We will explain the principle on an example of a general 3 ×3 kernel G3×3. Fig. 5a shows the topology for this type of kernel. Each element of the structure (Fig. 5b) contains two SISO modules [3]. SISO modules (in red) perform the soft inversion of the convolution itself. It is followed by the soft inversion of broadcasters (in black). The broadcasters are in the reference channel model responsible for input data branching and preparation of individual convolution regions. Data branching is described in a combination table Tab. 3, which controls the soft inversions of broadcasters in the decoding network. i7→A(B) IO ic[0]7→Ad··· ic[8]id7→Ad 0 0 ··· 0 0 1 1 ··· 1 1 Tab. 3. Table of combinations controlling the soft inversion of binary broadcaster. / k+ 1 k k−1 l−1l+ 1l 2  1  (a) The IDN topology with marked activation schedule. / SO(ˇc)[k−1, l+1,5] SI(ˇc)[k−1, l+1,5] SO(ˇc)[k−1, l, 6] SI(ˇc)[k−1, l, 6] SO(ˇc)[k, l, 2] SI(ˇc)[k, l, 2] SI(ˇc)[k, l, 3] SO(ˇc)[k, l, 3] SO(ˇc)[k, l, 0] SI(ˇc)[k, l, 0] SO(ˇc)[k, l, 8] SI(ˇc)[k, l, 8] SO(ˇc)[k, l, 1] SI(ˇc)[k, l, 1] SO(ˇc)[k+1, l+1,3] SI(ˇc)[k+1, l+1,3] SB(ˇq)[k, l] SF(ˇq)[k, l] SO(ˇc)[k, l+1,4] SI(ˇc)[k, l+1,4] SO(ˇc)[k+1, l−1,1] SI(ˇc)[k+1, l−1,1] SO(ˇc)[k, l, 6] SI(ˇc)[k, l, 6] SO(ˇc)[k+1, l, 2] SI(ˇc)[k+1, l, 2] SO(ˇc)[k, l, 7] SI(ˇc)[k, l, 7] SO(ˇc)[k, l−1,0] SI(ˇc)[k, l−1,0] SO(ˇc)[k, l, 8] SI(ˇc)[k, l, 8] SO(ˇc)[k, l, 5] SI(ˇc)[k, l, 5] SB(ˇ d)[k, l] SF(ˇ d)[k, l] SO(ˇc)[k−1, l−1,7] SI(ˇc)[k−1, l−1,7] SO(ˇc)[k, l, 4] SI(ˇc)[k, l, 4] B−1 f−1(.) (b) The IDN cell in the node [k,l]. / SO(ˇc)[k, l, 5] SI(ˇc)[k, l, 8] SI(ˇc)[k, l−1,0] SO(ˇc)[k, l, 4] SF(ˇ d)[k, l] SB(ˇ d)[k, l] SI(ˇc)[k−1, l−1,7] SI(ˇc)[k+1, l+1,3] SB(ˇq)[k, l] SF(ˇq)[k, l] SI(ˇc)[k, l+1,4] SO(ˇc)[k, l, 0] SO(ˇc)[k, l, 8] SO(ˇc)[k, l, 1] SO(ˇc)[k, l, 3] SO(ˇc)[k, l, 2] SI(ˇc)[k−1, l, 6] SI(ˇc)[k−1, l+1,5] SI(ˇc)[k+1, l−1,1] SO(ˇc)[k, l, 6] SI(ˇc)[k+1, l, 2] SO(ˇc)[k, l, 7] B f−1(.) (c) The simplified IDN cell in the node [k,l]. Fig. 5. Iterative decoding network marginalizing at the symbol level for a general ISI channel G3×3. 624 D. KEKRT, T. LUKE ˇ S, M. KL´ IMA, K. FLIEGEL, 2D ITERATIVE MAP DETECTION The activation schedule is parallel. It appears in purple in Fig. 5a. The activation takes place in all elements at once and consists of two steps. Firstly, metrics collected from the neighboring elements activate the soft inversion of the combinational logic. Secondly, soft inversions of broadcasters are performed followed by the distribution of output metrics into the neighboring elements. The whole process is repeated in the next iteration. The convolution is separated on the lowest possible level (symbol level). The network is flexible and can be well adjusted even for kernels with uncommon shapes. The disadvantage is lower performance caused by maximal decomposition of the convolution region and a large amount of connections. On the other hand, memory requirements are low because cardinality of the metrics of the state variables is the smallest possible. The variable topology can be implemented in a simplified version. It significantly reduces computational demands. A simplified element [11] is shown in Fig. 5c. The simplified implementation is based on the assumption that the center of the kernel has the dominant coefficient. The output of the reference ISI channel model is then mostly influenced by the center of the convolution region. Therefore, the element contains only one SISO module, which performs one marginalization. Only the state variable in the center is marginalized. The output is one metric and it is further distributed into the neighboring elements by a common broadcaster. Connections in the topology are the same, but compared to the layout in Fig. 5a the distribution of metrics is performed only in one direction (away from broadcaster). In the simplified version, only one marginalization is carried out (instead of nine). It hugely reduces computational complexity. Additionally, the soft inversion of the broadcaster is not necessary. The amount of connections is half that because the flow of information is only one-way. Memory requirements are also lower, but the performance deteriorates and the network can not be used to solve a more complex task (as for instance de-mosaicing, which is explained in the next chapter). 4.4 Layered Topologies The layered topology represents a combination of the two previously mentioned topologies. It can be used only if it is possible to decompose the convolution kernel into two orthogonal directions (as mentioned in 2.2). If the twodimensional convolution can be replaced by two consecutive one-dimensional convolutions, the reference channel model will have two layers. So as the iterative network has two layers. Each of these layers can be implemented using fixed or variable topology. The utilization of fixed topology leads to a detector where the lower layer is concatenated horizontally and the upper layer vertically (Fig. 6a). The SISO modules of the upper and lower level are shown in Fig. 6b and Fig. 6c, respectively. One iteration of the system is given by activation of the upper layer in a column by column manner. It is followed / k+ 1 k k−1 ll −1l+ 1 f−1 H(.) f−1 H(.) f−1 H(.) f−1 H(.) f−1 H(.) f−1 H(.)f−1 H(.)f−1 H(.) f−1 H(.) 2  1  f−1 V(.) f−1 V(.) f−1 V(.) f−1 V(.) f−1 V(.) f−1 V(.) f−1 V(.) f−1 V(.) f−1 V(.) (a) The IDN topology with marked activation schedule. / SF(ˇ d)[k, l] SB(ˇ d)[k, l] SO(ˇ R)[k, l] SI(ˇ R)[k, l] SF(ˇc)[k, l] SB(ˇc)[k, l] SI(ˇ R)[k, l + 1] SO(ˇ R)[k, l + 1] f−1 H(.) (b) The IDN cell in the node [k,l]on the bottom layer. / SI(ˇ C)[k+ 1, l] SO(ˇ C)[k+ 1, l]SI(ˇ C)[k, l] SO(ˇ C)[k, l] SF(ˇc)[k, l] SB(ˇc)[k, l] SF(ˇq)[k, l] SB(ˇq)[k, l] f−1 V(.) (c) The IDN cell in the node [k,l]on the top layer. Fig. 6. Layered iterative decoding network marginalizing at the symbol block level. by parallel activation of the lower level, which is done in a row by row manner. During the iteration process, the state metrics SI|O(ˇ R)and SI|O(ˇ C)on the individual layers and the “inter-layers” metrics SF|B(ˇc)are precised. If the variable topology is used, we obtain the structure in Fig. 7 (for the deconvolution of a 3 ×3 channel). The iteration process is then analogical. The activation runs from the upper layer and it is followed by the activation of the lower layer. The biggest advantage of the layered detector is the reduction of computational demands. Let us assume a decomposed kernel G(HV) 3×3. In the case of the standard (nonsimplified) variable topology and fixed topology, the cardinality of the SISO module combination table is card(AIO) = 29=512. For the layered topology, the cardinality of the lower layer is card(A(H) IO ) = 23=8 and on the upper layer card(A(V) IO ) = 63=216, which is 224 in total. The computational demand is decreased but as a consequence the performance is degraded. Basically, the performance decreases RADIOENGINEERING, VOL. 23, NO. 2, JUNE 2014 625 / k−1 k k+ 1 l−1l l + 1 3  4  2  1  Fig. 7. Layered iterative decoding network marginalizing at the symbol level for a general decomposition-able ISI channel G(HV) 3×3. with the increasing number of marginalization operations. The layered topology contains one extra marginalization that produces “inter-layers” metrics. 5. Application in Image Restoration The performance and capabilities of iterative detection will be demonstrated on a restoration of binary images obtained from a camera with CCD sensor. Thresholding that leads to a binary image is often an important first step of image analysis [19]. In many applications, images are in essence binary. They are blurred by the optical system and some noise is added during the process of image acquisition. The thresholding then can be seen as a reconstruction of the ideal image. The iterative decoding network reconstructs the desired image due to its deconvolution and noise suppression properties. The algorithm decides for each pixel whether its value will be black or white by maximizing the corresponding a posteriori probability that the pixel contains the useful signal (object in the foreground) or not (background noise etc.). We also demonstrate that if the image is undersampled, the iterative decoding network is able to calculate the missing information. It promises interesting applications in image interpolation. 5.1 System Model The ISI channel in this application is the objective of the camera. The diffraction of light prevents exact convergence of the light rays to a single point at the image plane. A sharp point on the object is blurred into a finite-sized spot in the image which is described by the point spread function (PSF) of the objective. The diffraction limited PSF of an ideal lens with a circular aperture can be described by the Airy disk [21]. The PSF of a real lens of a common camera can be well approximated by Gaussian model [22], PSF∆(x,y) = ∆ πe−∆(x2+y2)(13) where ∆determines the lobe width. Model does not consider the loss of light in the lens. Therefore ´∞ −∞´∞ −∞PSF∆(x,y)dxdy =1. Cells of the imaging sensor are assumed to be square-shaped, with no gaps between them and with a normalized length equal to 1. The amount of light that impacts a cell of the sensor is given by integration of the PSF over a particular square-shaped area. We obtain the Gaussian kernel {G(Gauss) L×L,∆}k,l=ˆl+1 lˆk+1 k PSF∆x−1 2,y−1 2dxdy (14) where Lis the size of the kernel, which determines a number of neighboring pixels affecting each other by blurring. Between the suppression of the main beam and lobe width of the PSF, there is a relation {G(Gauss) L×L,∆}0,0=erf(√∆/2)2. The kernel size Land the lobe width ∆reflect the properties of the camera objective. The camera sensor introduces random additive noise to the signal. The basic camera model takes into account three noise sources: thermal noise WT, readout noise WR and quantization noise WCof the A/D converter [14]. The signal at the output of the imaging sensor can be expressed as RC=N(ET) eQ+WT+WR ∆C (15) =S+WT+WR ∆C where ∆C=j2−NBN(FWC) ek(16) denotes the number of electrons per quantization step of the converter and N(ET) eis the number of electrons generated in the potential well at the maximum irradiation of the sensor (Q=1) during the exposure time TE. We expect an exposure time that does not cause the saturation of the sensor so that N(ET) e<N(FWC) ewhere N(FWC) eis the full well capacity (FWC) of the sensor. The signal is quantized at the output of the sensor, R=RC+WC(17) =dRCc using rounding d.cto the nearest integer value. Gateway of detector processes the signal Rusing the knowledge of statistical properties of the noise. Noise distributions are known. Under standard conditions, the Poisson distribution of the thermal noise can be approximated by normal distribution with the mean value µTand standard deviation √µT. The approximation holds true for values µT=40 [14]. The combination of these two noise sources leads to the Gaussian additive noise with the probability density pw(ξ,µw,σw) = 1 √2πσ exp(−(ξ−µw)2/2σ2 w). We integrate this continuous distribution over different quantization 626 D. KEKRT, T. LUKE ˇ S, M. KL´ IMA, K. FLIEGEL, 2D ITERATIVE MAP DETECTION steps from 0 up to the full well capacity. So we get discrete distribution, Pr(Cut) w(µ,σ)[n,N] =          1 2erfc2µ−1 √8σ,n=0 1 2erf2(n−µ)+1 √8σ−erf2(n−µ)−1 √8σ,0<n<N 1 2erfc2(N−µ)−1 √8σ,n=N (18) where Nwill be substituted for the highest quantization level. Thus, we obtain the transformation function of the gateway, {PF(q(i))[k,l]}i=(Pr(Cut) w ˆµT+N(ET) eq(i) ∆C , √ˆµT+ˆ σR ∆C[r[k,l],2NB−1]i (19) where ˆµTis an estimate of the mean value of the thermal noise, ˆ σRis an estimate of the standard deviation of the readout noise, NBis the number of bits of the A/D converter and q(i)are individual values from the alphabet Aq. 5.2 Examples of Dichromatic Image Restoration with Perfect CSI Knowledge Image binarization will be demonstrated using two types of convolution kernels. The first one is a 3 ×3 kernel with the transmission of the main ray 0.3. G(Gauss) 3×3,1.13 =   0.046 0.117 0.046 0.117 0.3000 0.117 0.046 0.117 0.046    (20) In the second case, a 5 ×5 reduced kernel is utilized with suppressed side beams and the transmission of the main ray 0.2, G(Gauss,R) 5×5,0.705 =            0.016 0.057 0.107 0.057 0.016 0.107 0.2000 0.107 0.016 0.057 0.107 0.057 0.016            .(21) Let us assume that the properties of the optical system on the detection part are known and that the detector works with perfect channel state information (CSI) about the values of the kernel. In this simulation, we model the image acquisition by the sensor iXon3 885 (Andor Technology). Real properties of this sensor were used for the following parameters N(FWC) e=30 ×103[e],∆C=14 [e],NB=11. For the chosen value N(ET) e=29×103[e], the sensor is excited up to 96.5 %. We leave ˆ σRas an independent variable, which determines the standard deviation ˆ σ=√ˆµT+ˆ σRof the resulting Gaussian noise. In the simulations, we used ˆµT=60 [e]. The transfer function of the gateway corresponds to the properties of the sensor in combination with the first model of the convolution kernel (20) (Fig. 8). For the second model of the convolution kernel, the transfer function is similar with typical “s-shaped deflection”. 010 20 30 40 0 500 1000 1500 2000 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 i CH: Histograms r[k, l] PF(q(i))[k,l] (a) ˆ σ=103[e]. 010 20 30 40 0 500 1000 1500 2000 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 i CH: Histograms r[k, l] PF(q(i))[k,l] (b) ˆ σ=316 [e]. Fig. 8. The examples of transfer functions of network gateway for Gaussian kernel 3×3, ∆=1.13 and various ˆ σ. Fig. 9 shows a simulated image of a QR code at the output of the A/D converter of a camera. Two cases are the outputs of a monochromatic camera and two are outputs from a camera with the Bayer mask (green filter pattern). Fig. 10 and Fig. 11 represent the reconstructed (binarized) output images together with their error rates after one (I=1), three (I=3) and five (I=5) iterations. In the case of 3 ×3 kernel, the image was reconstructed by the decoding network with a marginalization at the pixel block level and SyD technique. We used marginalization at the pixel level and PgD technique for the 5 ×5 kernel. The results show that the error rate falls rapidly with an increasing iteration number. Already the 3rd iteration provides satisfactory results. If the kernels are large, the level of noise should be lower or the detection may start to fail. It is caused by a relatively small bit-depth of the A/D converter. The quantization step ∆Cis then too coarse and the deconvolution task becomes ambiguous due to the larger cardinality of the al-