Joint Denoising / Compression of Image Contours via Shape Prior and Context Tree

TL;DR

This paper introduces a joint denoising and compression method for image contours using a shape prior and context tree, improving rate-distortion performance over separate approaches in bandwidth-constrained scenarios.

Contribution

It presents a theoretically justified joint denoising and compression framework with a novel dynamic programming algorithm and a compact context representation, outperforming separate methods.

Findings

Joint approach outperforms separate denoising and compression schemes.

Proposed DP algorithm achieves optimal solutions efficiently.

Experimental results confirm improved rate-distortion performance.

Abstract

With the advent of depth sensing technologies, the extraction of object contours in images---a common and important pre-processing step for later higher-level computer vision tasks like object detection and human action recognition---has become easier. However, acquisition noise in captured depth images means that detected contours suffer from unavoidable errors. In this paper, we propose to jointly denoise and compress detected contours in an image for bandwidth-constrained transmission to a client, who can then carry out aforementioned application-specific tasks using the decoded contours as input. We first prove theoretically that in general a joint denoising / compression approach can outperform a separate two-stage approach that first denoises then encodes contours lossily. Adopting a joint approach, we first propose a burst error model that models typical errors encountered in an observed string y of directional edges. We then formulate a rate-constrained maximum a posteriori (MAP) problem that trades off the posterior probability p(x'|y) of an estimated string x' given y with its code rate R(x'). We design a dynamic programming (DP) algorithm that solves the posed problem optimally, and propose a compact context representation called total suffix tree (TST) that can reduce complexity of the algorithm dramatically. Experimental results show that our joint denoising / compression scheme outperformed a competing separate scheme in rate-distortion performance noticeably.

Figures

Tables

Table 1. TABLE I: Results of Model Validation

Sequence	Noise	-Log-likelihood		AIC		Relative
		i.i.d	Prop	i.i.d	Prop	likelihood
Tsukuba	10%	67	64	134	128	3.93e-02
	30%	141	134	283	267	3.72e-04
Dude	10%	198	178	397	356	1.65e-09
	30%	378	348	756	696	8.93e-14

Equations

$$\underset{\hat{\mathbf{x}}\in \mathcal{S}}{\max }P(\hat{\mathbf{x}}|\mathbf{y})\text{mbox}s.t.R(\hat{\mathbf{x}})\le R_{\max }$$

$${\mathbf{x}}_{MAP}=\arg \underset{\hat{\mathbf{x}}\in \mathcal{S}}{\max }P(\hat{\mathbf{x}}|\mathbf{y})$$

$$\underset{\hat{\mathbf{x}}\in \mathcal{S}}{\max }P(\hat{\mathbf{x}}|{\mathbf{x}}_{MAP})\text{mbox}s.t.R(\hat{\mathbf{x}})\le R_{\max }$$

$$P(\hat{\mathbf{x}}|{\mathbf{x}}_{MAP})=\sum _{\hat{\mathbf{y}}}P(\hat{\mathbf{x}}|{\mathbf{x}}_{MAP},\hat{\mathbf{y}})P(\hat{\mathbf{y}}|{\mathbf{x}}_{MAP}).$$

$$P(\hat{\mathbf{y}}|{\mathbf{x}}_{MAP})=\frac{P({\mathbf{x}}_{MAP}|\hat{\mathbf{y}})P(\hat{\mathbf{y}})}{P({\mathbf{x}}_{MAP})}.$$

$$P({\mathbf{x}}_{MAP}|\hat{\mathbf{y}})=\{\begin{array}{ll}1,& \text{if }{\mathbf{x}}_{MAP}=\arg \underset{\tilde{\mathbf{x}}\in \mathcal{S}}{\max }P(\tilde{\mathbf{x}}|\hat{\mathbf{y}})\\ 0,& \text{otherwise}\end{array}.$$

$$\begin{array}{rl}& \sum _{\hat{\mathbf{y}}}P(\hat{\mathbf{x}}|{\mathbf{x}}_{MAP},\hat{\mathbf{y}})P({\mathbf{x}}_{MAP}|\hat{\mathbf{y}})P(\hat{\mathbf{y}})\\ & =\sum _{\hat{\mathbf{y}}|{\mathbf{x}}_{MAP}=\arg \underset{\tilde{\mathbf{x}}\in \mathcal{S}}{\max }P(\tilde{\mathbf{x}}|\hat{\mathbf{y}})}P(\hat{\mathbf{x}}|{\mathbf{x}}_{MAP},\hat{\mathbf{y}})P(\hat{\mathbf{y}})\\ & =\sum _{\hat{\mathbf{y}}|{\mathbf{x}}_{MAP}=\arg \underset{\tilde{\mathbf{x}}\in \mathcal{S}}{\max }P(\tilde{\mathbf{x}}|\hat{\mathbf{y}})}P(\hat{\mathbf{x}}|\hat{\mathbf{y}})P(\hat{\mathbf{y}})\end{array}$$

$$\underset{\hat{\mathbf{x}}\in \mathcal{S}}{\max }\sum _{\hat{\mathbf{y}}|{\mathbf{x}}_{MAP}=\arg \underset{\tilde{\mathbf{x}}\in \mathcal{S}}{\max }P(\tilde{\mathbf{x}}|\hat{\mathbf{y}})}P(\hat{\mathbf{x}}|\hat{\mathbf{y}})P(\hat{\mathbf{y}})\text{mbox}s.t.R(\hat{\mathbf{x}})\le R_{\max }.$$

$$\underset{\hat{\mathbf{x}}\in \mathcal{S}}{\min }-\log P(\hat{\mathbf{x}}|\mathbf{y})+\lambda R(\hat{\mathbf{x}})$$

$$P(\mathbf{x}|\mathbf{y})=\frac{P(\mathbf{y}|\mathbf{x})P(\mathbf{x})}{P(\mathbf{y})}$$

$$\begin{array}{rl}& (1-p{)}^{l_0(0)}\prod _{k=1}^{K}p(1-q_1{)}^{l_1(k)-1}{q}_1(1-q_2{)}^{l_2(k)-1}{q}_2(1-p{)}^{l_0(k)-1}\end{array}$$

$$\begin{array}{rl}& -\log P(\mathbf{y}|\mathbf{x})=\\ & -K(\log p+\log q_1+\log q_2)-(\Gamma -K)\log (1-p)\\ & -(\Lambda -K)\log (1-q_1)-(\Delta -K)\log (1-q_2)\end{array}$$

$$\begin{array}{rl}& -\log P(\mathbf{y}|\mathbf{x})\\ & \approx -K(\log p+\log q_1+\log q_2)\\ & -(\Lambda -K)\log (1-q_1)-(\Delta -K)\log (1-q_2)\\ & =-K(\underset{-c_0}{\underbrace{\log p+\log \frac{q_1}{1-q_1}+\log \frac{q_2}{1-q_2}}})\\ & -\Lambda \underset{-c_1}{\underbrace{\log (1-q_1)}}-\Delta \underset{-c_2}{\underbrace{\log (1-q_2)}}\end{array}$$

$$\begin{array}{rl}& -\log P(\mathbf{y}|\mathbf{x})\approx (c_0+c_2)K+c_1\Lambda +c_2{\Delta }^{\prime }\end{array}$$

$$P(\mathbf{x})=\exp \left\{-\beta \sum _{i=D_s+1}^{L_{\mathbf{x}}}s({\mathbf{x}}_{i-D_s}^i)\right\}$$

$$\begin{array}{rl}& \underset{0\le k\le L_{\mathbf{w}}}{\max }\left\{\frac{|(m_k-m_0)(n_{L_{\mathbf{w}}}-n_0)-(n_k-n_0)(m_{L_{\mathbf{w}}}-m_0)|}{\sqrt{(m_{L_{\mathbf{w}}}-m_0{)}^2+(n_{L_{\mathbf{w}}}-n_0{)}^2}}\right\}\end{array}$$

$$\begin{array}{rl}-\log P(\hat{\mathbf{x}}|\mathbf{y})=& (c_0+c_2)K+c_1\Lambda +c_2{\Delta }^{\prime }\\ & +\beta \sum _{i=D_s+1}^{L_{\hat{\mathbf{x}}}}s({\hat{\mathbf{x}}}_{i-D_s}^i)\end{array}$$

$$R(\hat{\mathbf{x}})=-\sum _{i=1}^N{\log }_{2}P({\hat{x}}_i|{\hat{\mathbf{x}}}_1^{i-1})$$

$$N(\mathbf{u})=\sum _{m=1}^M\sum _{i=1}^{L_{\mathbf{x}(m)}-|\mathbf{u}|+1}1\left(\mathbf{x}(m{)}_i^{i+|\mathbf{u}|-1}=\mathbf{u}\right)$$

$$\begin{array}[]{cc}P(x|\mathbf{u})=\frac{N(x\mathbf{u})}{\delta+N(\mathbf{u})}\end{array}$$

$$P({\hat{x}}_i|{\hat{\mathbf{x}}}_1^{i-1},\text{context model})=P({\hat{x}}_i|\mathbf{w})$$

$$P({\hat{x}}_i|{\hat{\mathbf{x}}}_1^{i-1},\mathcal{T})=P({\hat{x}}_i|\mathbf{w}).$$

$$P({\hat{x}}_i|\mathbf{w})=\{\begin{array}{ll}\frac{N({\hat{x}}_i\mathbf{w})}{|{\mathcal{A}}_{\mathbf{w}}|+N(\mathbf{w})},& \text{if }{\hat{x}}_i\in {\mathcal{A}}_{\mathbf{w}}\\ \frac{|{\mathcal{A}}_{\mathbf{w}}|}{|{\mathcal{A}}_{\mathbf{w}}|+N(\mathbf{w})}\cdot P({\hat{x}}_i|{\mathbf{w}}_2^{|\mathbf{w}|}),& \text{otherwise}\end{array}.$$

$$\begin{array}{rl}J(\hat{\mathbf{x}})=& -\log P(\mathbf{y}|\hat{\mathbf{x}})-\beta \log P(\hat{\mathbf{x}})+\lambda R(\hat{\mathbf{x}})\\ \approx & (c_0+c_2)K+c_1\Lambda +c_2{\Delta }^{\prime }+\beta \sum _{i=D_s+1}^{L_{\hat{\mathbf{x}}}}s({\hat{\mathbf{x}}}_{i-D_s}^i)\\ & -\lambda \sum _{i=D+1}^{L_{\hat{\mathbf{x}}}}{\log }_{2}P({\hat{x}}_i|{\hat{\mathbf{x}}}_{i-D}^{i-1}).\end{array}$$

$$G_i({\hat{\mathbf{x}}}_{i-D}^{i-1},\mathbf{e},j-1)=\underset{{\hat{x}}_i\in \mathcal{A}}{\min }\{\begin{array}{ll}f({\hat{\mathbf{x}}}_{i-D}^i)+1(j<L_{\mathbf{y}})G_{i+1}({\hat{\mathbf{x}}}_{i-D+1}^i,v(\mathbf{e},{\hat{x}}_i),j),& \text{if}{\hat{x}}_i=y_j\\ (c_0+c_2)+f({\hat{\mathbf{x}}}_{i-D}^i)+B_{i+1}({\hat{\mathbf{x}}}_{i-D+1}^i,v(\mathbf{e},{\hat{x}}_i),j),& \text{otherwise}\end{array}$$

$$B_i({\hat{\mathbf{x}}}_{i-D}^{i-1},\mathbf{e},j-1)=\underset{{\hat{x}}_i\in \mathcal{A}}{\min }\{\begin{array}{ll}{c}_2(k-j)+f({\hat{\mathbf{x}}}_{i-D}^i)+G_{i+1}({\hat{\mathbf{x}}}_{i-D+1}^i,v(\mathbf{e},{\hat{x}}_i),k),& \text{if}\exists k,v(\mathbf{e},{\hat{x}}_i)={\mathbf{e}}_{\mathbf{y}}(k)\\ c_1+f({\hat{\mathbf{x}}}_{i-D}^i)+B_{i+1}({\hat{\mathbf{x}}}_{i-D+1}^i,v(\mathbf{e},{\hat{x}}_i),j),& \text{otherwise}\end{array}$$

$$f({\hat{\mathbf{x}}}_{i-D}^i)=\beta s({\hat{\mathbf{x}}}_{i-D_s}^i)-\lambda {\log }_{2}P({\hat{x}}_i|{\hat{\mathbf{x}}}_{i-D}^{i-1}).$$

$$D(\hat{\mathbf{x}},\mathbf{x})=\sum _{i=1}^{L_{\hat{\mathbf{x}}}}{d}^2({\mathbf{e}}_{\hat{\mathbf{x}}}(i),\mathbf{x}).$$

$$d({\mathbf{e}}_{\hat{\mathbf{x}}}(i),\mathbf{x})=\underset{1\le j\le L_{\mathbf{x}}}{\min }|{\hat{m}}_i-m_j|+|{\hat{n}}_i-n_j|$$

$$AIC=2k-2\ln \mathcal{L}$$

Full text

I Introduction

Advances in depth sensing technologies like Microsoft Kinect 2.0 means that depth images—per pixel distances between physical objects in a 3D scene and the camera—can now be captured easily and inexpensively. Depth imaging has in turn eased the extraction of object contours in a captured image, which was traditionally a challenging computer vision problem [1]. Detected contours can be used to facilitate recent advanced image coding schemes; if object contours are compressed efficiently as side information (SI), they can enable new edge-adaptive techniques such as graph Fourier transform (GFT) coding [2, 3] and motion prediction of arbitrarily shaped blocks [4]. Further, coded object contours can also be transmitted to a central cloud for computationally expensive application-specific tasks such as object detection or human action recognition [5], at a much lower coding cost than the original captured depth video.

Unfortunately, captured depth images are often corrupted by acquisition noise, and hence unavoidably the detected contours also contain errors. In this paper, we propose to jointly denoise and compress detected object contours in images. First, we prove theoretically that in general a joint denoising / compression approach outperforms a two-stage separate approach that first denoises an observed contour then compresses the denoised version lossily. Adopting a joint approach, we first propose a burst error model that captures unique characteristics of typical errors encountered in detected contours. We then formulate a rate-constrained maximum a posteriori (MAP) problem that trades off the posterior probability $P(\hat{\mathbf{x}}|\mathbf{y})$ of an estimated contour $\hat{\mathbf{x}}$ given observed contour $\mathbf{y}$ with its code rate $R(\hat{\mathbf{x}})$. Given our burst error model, we show that the negative log of the likelihood $P(\mathbf{y}|\mathbf{x})$ can be written intuitively as a simple sum of burst error events, error symbols and burst lengths. Further, we construct a geometric prior $P(\mathbf{x})$ stating intuitively that contours are more likely straight than curvy.

We design a dynamic programming (DP) [6] algorithm that solves the posed problem optimally, and propose a compact context representation called total suffix tree (TST) that can reduce the algorithm complexity dramatically. Experimental results show that our joint denoising / compression scheme outperformed a competing separate scheme in RD performance noticeably. We note that, to the best of our knowledge, we are the first in the literature111An early version of this work appeared as a conference paper [7]. to formally address the problem of joint denoising / compression of detected image contours.

The outline of the paper is as follows. We first overview related works in Section II. We pose our rate-constrained MAP problem in Section III, and define the corresponding error and rate terms in Section IV and V, respectively. We describe our proposed optimal algorithm in Section VI. Experimental results are presented in Section VII, and we conclude in Section VIII.

II Related Work

II-A Contour Coding

II-A1 Lossless Contour Coding

Most works in lossless contour coding [4, 8, 9, 10, 11, 12, 13, 14, 15] first convert an image contour into a chain code [16]: a sequence of symbols each representing one of four or eight possible absolute directions on the pixel grid. Alternatively, a differential chain code (DCC) [8] that specifies relative directions instead can be used. DCC symbols are entropy-coded using either Huffman [9] or arithmetic coding [10] given symbol probabilities. The challenge is to estimate conditional probabilities for DCC symbols given a set of training data; this is the main problem we discuss in Section V.

[17, 4] propose a linear geometric model to estimate conditional probabilities of the next DCC symbol. In summary, given a window of previous edges, a line-of-best-fit that minimizes the sum of distances to the edges’ endpoints is first constructed. Then the probability of a candidate direction for the next symbol is assumed inversely proportional to the angle difference between the direction and the fitted line. This scheme is inferior in estimating symbol probabilities compared to context models, because there are only a few possible angle differences for a small number of previous edges, limiting the expressiveness of the model.

An alternative approach is context modeling: given a window of $l$ previous symbols (context) $\mathbf{x}^{i-1}_{i-l}$, compute the conditional probability $P(x_{i}|\mathbf{x}^{i-1}_{i-l})$ of the next symbol $x_{i}$ by counting the number of occurrences of $\mathbf{x}^{i-1}_{i-l}$ followed by $x_{i}$ in the training data. In [18, 10, 11, 12, 13], Markov models of fixed order up to eight are used for lossless coding. However, in applications where the training data is limited, there may not be enough occurrences of $\mathbf{x}^{i-1}_{i-l}$ to reliably estimate the conditional probabilities.

Variable-length Context Tree (VCT) [19, 20] provides a more flexible approach for Markov context modeling by allowing each context to be of variable length. There are many ways to construct a VCT; e.g., Lempel-Ziv-78 (LZ78) [21] and prediction by partial matching (PPM) [22, 23]. LZ78 constructs a dictionary from scratch using the input data directly as training data. The probability estimation quality varies depending on the order of first appearing input symbols. PPM considers all contexts restricted by a maximum length with non-zero occurrences in the training data when building a VCT. PPM has an efficient way to deal with the zero frequency problem [20], where the input symbol along with the context does not occur in the training data, by reducing the context length. PPM is widely used for lossless contour coding [14, 15] because of its efficiency. In this paper, we use PPM for contour coding, as described in Section V.

II-A2 Lossy Contour Coding

Lossy contour coding approaches can be classified into two categories: DCC-based [24, 25] and vertex-based [26, 27, 28] approaches. In [24, 25], the DCC strings are approximated by using line smoothing techniques. In contrast, vertex-based approaches select representative key points along the contour for coding at the encoder and interpolation at the decoder. Because vertex-based approaches are not suitable for lossless contour coding and we consider joint denoising / compression of contours for a wide range of bitrates, we choose a DCC-based approach and compute the optimal rate-constrained MAP contour estimate for coding using PPM.

II-B Contour Denoising

Contour noise removal is considered in [29, 4, 30, 31, 32, 33]. In [29], the noisy pixels along the chain code are first removed, then the processed chain code is decomposed into a set of connected straight lines. In [4], a chain code is simplified by removing “irregularities”, which are predefined non-smooth edge patterns. Both [29] and [4] are contour smoothing techniques, which does not guarantee the similarity between the original contour and the smoothed contour. Sparse representation and non-linear filter are used in [30] and [31] respectively to enhance the contours in a binary document image, with the objective of removing the noise in the binary image. [32, 33] denoise the image contours using Gaussian filter by averaging the coordinates of the points on the contours. Gaussian filter is widely used for removing noise from signal, but the Gaussian model does not specifically reflect the statistics of the noise on the image contours. In contrast, we propose a burst error model to describe the errors on the image contours, as discussed in Section IV.

Considering the problem of both denoising and compression of the image contours, the denoised contour by [29, 4, 30, 31, 32, 33] may require a large encoding overhead. We perform a MAP estimate to denoise an observed contour subject to a rate constraint. We will show in Section VII that we outperform a separate denoising / compression approach.

III Problem Formulation

We first describe our contour denoising / compression problem at a high level. We assume that one or more object contours in a noise-corrupted image have first been detected, for example, using a method like gradient-based edge detection [4]. Each contour is defined by an initial point and a following sequence of connected “between-pixel” edges on a 2D grid that divide pixels in a local neighbourhood into two sides. As an example, three contours in one video frame of Dude are drawn in Fig. 1. For each detected (and noise-corrupted) contour, the problem is to estimate a contour that is highly probable and requires few bits for encoding.

We discuss two general approaches for this problem—joint approach and separate approach—and compare and contrast the two. Adopting the joint approach, we formulate a rate-constrained MAP problem.

III-A Joint and Separate Approaches

III-A1 Definitions

We first convert each contour into a differential chain code (DCC)[8]. A length-$L_{\mathbf{z}}$ contour $\mathbf{z}$ can be compactly described as a symbol string denoted by $[z_{1},z_{2},\cdots,z_{L_{\mathbf{z}}}]$. As shown in Fig. 2, the first symbol $z_{1}$ is chosen from a size-four alphabet, $\mathcal{D}=\{\texttt{N},\texttt{E},\texttt{S},\texttt{W}\}$, specifying the four absolute directions north, east, south and west with respect to a 2D-grid start point, i.e., $z_{1}\in\mathcal{D}$. Each subsequent DCC symbol $z_{i}$, $i\geq 2$, is chosen from a size-three alphabet, $\mathcal{A}=\{\texttt{l},\texttt{s},\texttt{r}\}$, specifying the three relative directions left, straight and right with respect to the previous symbols $z_{i-1}$, i.e., $z_{i}\in\mathcal{A}$ for $i\geq 2$. Denote by $\mathbf{z}^{j}_{i}=[{z}_{j},{z}_{j-1},\ldots,{z}_{i}]$, $i<j$ and $i,j\in\mathbb{Z}^{+}$, a sub-string of length $j-i+1$ from the $i$-th symbol $z_{i}$ to the $j$-th symbol $z_{j}$ in reverse order.

Alternatively, one can represent the $i$-th edge of contour $\mathbf{z}$ geometrically on the 2D grid as $\mathbf{e}_{\mathbf{z}}(i)=\{(m_{i},n_{i}),d_{i}\}$, where $(m_{i},n_{i})$ is the 2D coordinate of the ending point of the edge and $d_{i}\in\mathcal{D}$ is the absolute direction, as shown in Fig. 2(a). This representation is used later when describing our optimization algorithm.

III-A2 Joint Problem Formulation

Denote by $\mathbf{y}\in\mathcal{S}$ and $\mathbf{x}\in\mathcal{S}$ the observed and ground truth DCC strings respectively, where $\mathcal{S}$ is the space of all DCC strings of finite length. As done in [34, 35] for joint denoising / compression of multiview depth images222We note that though [34, 35] also formulated a rate-constrained MAP problem for the joint denoising / compresion problem, there was no theoretical analysis on why a joint approach is better in general., we follow a rate-constrained MAP formulation and define the objective as finding a string $\hat{\mathbf{x}}\in\mathcal{S}$ that maximizes the posterior probability $P(\hat{\mathbf{x}}|\mathbf{y})$, subject to a rate constraint on chosen $\hat{\mathbf{x}}$:

[TABLE]

where $R(\hat{\mathbf{x}})$ is the bit count required to encode string $\hat{\mathbf{x}}$, and $R_{\max}$ is the bit budget.

(1) essentially seeks one solution to both the estimation problem and the compression problem at the same time. The estimated rate-constrained DCC string from (1) is then losslessly compressed at the encoder.

III-A3 Separate Problem Formulation

An alternative approach is to first solve an unconstrained MAP estimation sub-problem, then lossily compress the estimated DCC string subject to the rate constraint as follows:

[TABLE]

[TABLE]

where $\mathbf{x}_{MAP}$ is an optimal solution to the unconstrained estimation sub-problem (2). The lossy compression sub-problem (3) aims to find a DCC string $\mathbf{\hat{x}}$ which is closest to $\mathbf{x}_{MAP}$ in probability given bit budget $R_{\max}$. In the sequel we call (1) the joint approach, and (2) and (3) the separate approach.

Since the contour denoising / compression problem is inherently a joint estimation / compression problem, the solution to the joint approach should not be worse than the solution to the separate approach. The question is under what condition(s) the joint approach is strictly better than the separate approach.

III-B Comparison between Joint and Separate Approaches

III-B1 Algebraic Comparison

To understand the condition(s) under which the joint approach is strictly better than the separate approach, we first rewrite $P(\hat{\mathbf{x}}|\mathbf{x}_{MAP})$ in (3) using the law of total probability by introducing an observation $\hat{\mathbf{y}}$:

[TABLE]

$P(\hat{\mathbf{y}}|\mathbf{x}_{MAP})$ can be further rewritten using Bayes’ rule:

[TABLE]

From (2), $P(\mathbf{x}_{MAP}|\hat{\mathbf{y}})=1$ if $\mathbf{x}_{MAP}$ is the MAP solution given this observation $\hat{\mathbf{y}}$; otherwise, $P(\mathbf{x}_{MAP}|\hat{\mathbf{y}})=0$:

[TABLE]

Combining (2) and (4) to (6), $P(\hat{\mathbf{x}}|\mathbf{x}_{MAP})$ is rewritten as follows:

[TABLE]

where $P(\hat{\mathbf{x}}|\mathbf{x}_{MAP},\hat{\mathbf{y}})=P(\hat{\mathbf{x}}|\hat{\mathbf{y}})$, since $\mathbf{x}_{MAP}$ is a function of $\hat{\mathbf{y}}$, i.e., $\mathbf{x}_{MAP}=\arg\,\underset{\tilde{\mathbf{x}}\in\mathcal{S}}{\max}\ P(\tilde{\mathbf{x}}|\hat{\mathbf{y}})$. We ignore $P(\mathbf{x}_{MAP})$ in (5), because $\mathbf{x}_{MAP}$ (and hence $P(\mathbf{x}_{MAP})$) is fixed and does not affect the maximization.

Thus, we can rewrite the separate approach as:

[TABLE]

Comparing (8) to the joint approach (1), we can conclude the following: if there is only one observation $\hat{\mathbf{y}}$ that corresponds to $\mathbf{x}_{MAP}$, i.e., $\mathbf{x}_{MAP}=\arg\,\underset{\tilde{\mathbf{x}}\in\mathcal{S}}{\max}\ P(\tilde{\mathbf{x}}|\hat{\mathbf{y}})$, then the solutions to the separate and joint approaches are the same—they both maximize $P(\hat{\mathbf{x}}|\mathbf{y})$ subject to $R(\hat{\mathbf{x}})\leq R_{\max}$. This makes intuitive sense: if there is a one-to-one correspondence between $\mathbf{x}_{MAP}$ and observation $\mathbf{y}$, then knowing $\mathbf{x}_{MAP}$ is as good as knowing $\mathbf{y}$. On the other hand, if there are more than one observation $\hat{\mathbf{y}}$ corresponding to $\mathbf{x}_{MAP}$, i.e., there are more than one term in the sum in (8), then the solutions between the two approaches can be far apart.

III-B2 Geometric Comparison

To develop intuition on the difference between the joint and separate approaches, we examine the problem from a geometric viewpoint by representing each DCC string as a point in a high-dimensional space. By Bayes’ rule, maximization of the posterior $P(\hat{\mathbf{x}}|\mathbf{y})$ is the same as maximization of the product of the likelihood $P(\mathbf{y}|\hat{\mathbf{x}})$ and the prior $P(\hat{\mathbf{x}})$, or minimization of $-\log P(\mathbf{y}|\hat{\mathbf{x}})-\log P(\hat{\mathbf{x}})$. As shown in Fig. 3, each point in the space denotes a DCC string333We can map a DCC string into a point in a $N$-dimensional space as follows. First divide the string on the 2D grid into $(N-1)/2$ segments of equal length. The $x$- and $y$-coordinates of the end point of each segment on the 2D grid becomes an index for each segment. The resulted $N-1$ index numbers, plus a number to denote the length of each segment, becomes an $N$-dimension coordinate in an $N$-dimensional space.. Specifically, $\mathbf{y}$ and $\mathbf{y}^{\prime}$ in Fig. 3(a) denote two different observations that would yield the same unconstrained MAP solution $\mathbf{x}_{MAP}$, i.e., $\mathbf{x}_{MAP}=\arg\max_{\hat{\mathbf{x}}}P(\hat{\mathbf{x}}|\mathbf{y})=\arg\max_{\hat{\mathbf{x}}}P(\hat{\mathbf{x}}|\mathbf{y}^{\prime})$.

Denote by $\mathbf{x}_{0}$ the DCC string with the largest prior $P(\mathbf{x}_{0})$, e.g., the string corresponding to the “straightest” contour. For illustration, we draw isolines centered at $\mathbf{y}$ (or $\mathbf{y}^{\prime}$), each denoting the set of strings that have the same metric distance $D_{l}(\mathbf{y},\mathbf{x})$—negative log likelihood $-\log P(\mathbf{y}|\hat{\mathbf{x}})$—given observation $\mathbf{y}$ (or $\mathbf{y}^{\prime}$). Similarly, we draw isolines centered at $\mathbf{x}_{0}$ to denote strings with the same metric distance $D_{p}(\mathbf{x}_{0},\mathbf{x})$—negative log prior $-\log P(\hat{\mathbf{x}})$. Thus, maximizing $P(\hat{\mathbf{x}}|\mathbf{y})$ is equivalent to identifying a point $\mathbf{x}$ in the $N$-dimensional space that minimizes the sum of distance $D_{l}(\mathbf{y},\mathbf{x})$ and $D_{p}(\mathbf{x}_{0},\mathbf{x})$.

The joint approach is then equivalent to finding a point $\hat{\mathbf{x}}$ inside the isoline defined by the rate constraint $R(\mathbf{\hat{x}})=R_{\max}$—feasible rate region—to minimize the sum of the two aforementioned distances. On the other hand, the separate approach is equivalent to first finding the point $\mathbf{x}_{MAP}$ to minimize $D_{l}(\mathbf{y},\mathbf{x}_{MAP})+D_{p}(\mathbf{x}_{0},\mathbf{x}_{MAP})$, then find a point $\hat{\mathbf{x}}$ inside the feasible rate region to minimize $D_{l}(\mathbf{x}_{MAP},\hat{\mathbf{x}})$, i.e., the distance from $\hat{\mathbf{x}}$ to $\mathbf{x}_{MAP}$.

In Fig. 3(a), $\mathbf{x}_{c}^{\prime}$ is the solution to the separate approach, which is the point bounded by the rate constraint with minimum distance $D_{l}(\mathbf{x}_{MAP},\mathbf{x}_{c}^{\prime})$ to $\mathbf{x}_{MAP}$. Note that the solution to the joint approach denoted by $\mathbf{x}_{c}$ may be different from $\mathbf{x}_{c}^{\prime}$; i.e., the sum of the two distances $D_{l}(\mathbf{y},\mathbf{x}_{c})+D_{p}(\mathbf{x}_{0},\mathbf{x}_{c})$ from $\mathbf{x}_{c}$ is smaller than the distance sum from $\mathbf{x}_{c}^{\prime}$.

Geometrically, if the major axes444Major axis of an ellipse is its longer diameter. of the isolines centered at $\mathbf{y}$ and $\mathbf{x}_{0}$ are aligned in a straight line between $\mathbf{y}$ and $\mathbf{x}_{0}$ as shown in Fig. 3(b), then: i) $\mathbf{x}_{MAP}$ must reside on the line, and ii) there is only one corresponding observation $\mathbf{y}$ for a given $\mathbf{x}_{MAP}$. The first claim is true because any point not on the straight line between $\mathbf{y}$ and $\mathbf{x}_{0}$ will have a distance sum strictly larger than a point on the straight line and thus cannot be the optimal $\mathbf{x}_{MAP}$. Given the first claim, the second claim is true because $\mathbf{y}$ must be located at a unique point on an extrapolated line from $\mathbf{x}_{0}$ to $\mathbf{x}_{MAP}$ such that $\mathbf{x}_{MAP}=\arg\min_{\mathbf{x}}D_{l}(\mathbf{y},\mathbf{x})+D_{p}(\mathbf{x}_{0},\mathbf{x})$. In this special case, the solutions to the joint and separate approaches are the same. Clearly, in the general case (as shown in Fig. 3(a)) the axes of the isolines are not aligned, and thus in general the joint approach would lead to a better solution.

III-C Rate-Constrained Maximum a Posteriori (MAP) Problem

Since the condition of one-to-one correspondence between $\mathbf{y}$ and $\mathbf{x}_{MAP}$ is not satisfied in general, we formulate our rate-constrained MAP problem based on the joint approach. Instead of (1), one can solve the corresponding Lagrangian relaxed version instead [36]:

[TABLE]

where Lagrangian multiplier $\lambda$ is chosen so that the optimal solution $\hat{\mathbf{x}}$ to (9) has rate $R(\hat{\mathbf{x}})\leq R_{\max}$. [36, 37] discussed how to select an appropriate $\lambda$. We focus on how (9) is solved for a given $\lambda>0$.

We call the first term and the second term in (9) error term and rate term, respectively. The error term measures, in some sense, the distance between the estimated DCC string $\hat{\mathbf{x}}$ and the observed DCC string $\mathbf{y}$. Note that this error term is not the distance between the estimated DCC string $\hat{\mathbf{x}}$ and the ground truth DCC string $\mathbf{x}$, since we have no access to ground truth $\mathbf{x}$ in practice. Thus we seek an estimate $\hat{\mathbf{x}}$ which is somehow close to the observation $\hat{\mathbf{y}}$, while satisfying some priors of the ground truth signal. The estimate $\hat{\mathbf{x}}$ is also required to have small rate term $R(\hat{\mathbf{x}})$.

Next, we first specify the error term by using a burst error model and a geometric prior in Section IV. Then we will specify the rate term using a general variable-length context tree (VCT) model [20] in Section V.

IV Error Term

We first rewrite the posterior $P(\mathbf{x}|\mathbf{y})$ in (9) using Bayes’ Rule:

[TABLE]

where $P(\mathbf{y}|\mathbf{x})$ is the likelihood of observing DCC string $\mathbf{y}$ given ground truth $\mathbf{x}$, and $P(\mathbf{x})$ is the prior which describes a priori knowledge about the target DCC string. We next describe an error model for DCC strings, then define likelihood $P(\mathbf{y}|\mathbf{x})$ and prior $P(\mathbf{x})$ in turn.

IV-A Error Model for DCC String

Assuming that pixels in an image are corrupted by a small amount of independent and identically distributed (iid) noise, a detected contour will occasionally be shifted from the true contour by one or two pixels. However, the computed DCC string from the detected contour will experience a sequence of wrong symbols—a burst error. This is illustrated in Fig. 4(a), where the left single erred pixel (in red) resulted in two erred symbols in the DCC string. The right single error pixel also resulted in a burst error in the observed string, which is longer than the original string. Based on these observations, we propose our DCC string error model as follows.

We define a three-state Markov model as illustrated in Fig. 4(b) to model the probability of observing DCC string $\mathbf{y}$ given original string $\mathbf{x}$. State 0 is the good state, and burst error state 1 and burst length state 2 are the bad states. $p$, $q_{1}$ and $q_{2}$ are the transition probabilities from state 0 to 1, 1 to 2, and 2 to 0, respectively. Note that state 1 cannot transition directly to 0, and likewise state 2 to 1 and 0 to 2.

Starting at good state 0, each journey to state 1 then to 2 then back to 0 is called a burst error event. From state 0, each self-loop back to 0 with probability $1-p$ means that the next observed symbol $y_{i}$ is the same as $x_{i}$ in original $\mathbf{x}$. A transition to burst error state 1 with probability $p$, and each subsequent self-loop with probability $1-q_{1}$, mean observed $y_{i}$ is now different from $x_{i}$. A transition to burst length state 2 then models the length increase in observed $\mathbf{y}$ over original $\mathbf{x}$ due to this burst error event: the number of self-loops taken back to state 2 is the increase in number of symbols. A return to good state 0 signals the end of this burst error event.

IV-B Likelihood Term

Given the three-state Markov model, we can compute likelihood $P(\mathbf{y}|\mathbf{x})$ as follows. For simplicity, we assume that $\mathbf{y}$ starts and ends at good state 0. Denote by $K$ the total number of burst error events in $\mathbf{y}$ given $\mathbf{x}$. Further, denote by $l_{1}(k)$ and $l_{2}(k)$ the number of visits to state 1 and 2 respectively during the $k$-th burst error event. Similarly, denote by $l_{0}(k)$ the number of visits to state 0 after the $k$-th burst error event. We can then write the likelihood $P(\mathbf{y}|\mathbf{x})$ as:

[TABLE]

For convenience, we define the total number of visits to state 0, 1 and 2 as $\Gamma=\sum_{k=0}^{K}l_{0}(k)$, $\Lambda=\sum_{k=1}^{K}l_{1}(k)$ and $\Delta=\sum_{k=1}^{K}l_{2}(k)$, respectively. We can then write the negative log of the likelihood as:

[TABLE]

Assuming that burst errors are rare events, $p$ is small and $\log(1-p)\approx 0$. Hence:

[TABLE]

Thus $-\log P(\mathbf{y}|\mathbf{x})$ simplifies to:

[TABLE]

where $\Delta^{\prime}=\Delta-K$ is the length increase in observed $\mathbf{y}$ compared to $\mathbf{x}$ due to the $K$ burst error events555Length increase of observed $\mathbf{y}$ due to $k$-th burst error event is $l_{2}(k)-1$.. (14) states that the negative log of the likelihood is a linear sum of three terms: i) the number of burst error events $K$; ii) the number of error corrupted symbols $\Lambda$; and iii) the length increase $\Delta^{\prime}$ in observed string $\mathbf{y}$. This agrees with our intuition that more error events, more errors and more deviation in DCC length will result in a larger objective value in (9). We will validate this error model empirically in our experiments in Section VII.

IV-C Prior Term

Similarly to [4, 38], we propose a geometric shape prior based on the assumption that contours in natural images tend to be more straight than curvy. Specifically, we write prior $P(\mathbf{x})$ as:

[TABLE]

where $\beta$ and $D_{s}$ are parameters. $s(\mathbf{x}_{i-D_{s}}^{i})$ measures the straightness of DCC sub-string $\mathbf{x}_{i-D_{s}}^{i}$. Let $\mathbf{w}$ be a DCC string of length $D_{s}+1$, i.e., $L_{\mathbf{w}}=D_{s}+1$. Then $s(\mathbf{w})$ is defined as the maximum Euclidean distance between any coordinates of edge $\mathbf{e}_{\mathbf{w}}(k),0\leq k\leq L_{\mathbf{w}}$ and the line connecting the first point $(m_{0},n_{0})$ and the last point $(m_{L_{\mathbf{w}}},n_{L_{\mathbf{w}}})$ of $\mathbf{w}$ on the 2D grid. We can write $s(\mathbf{w})$ as:

[TABLE]

Since we compute the sum of a straightness measure for overlapped sub-strings, the length of each sub-string $\mathbf{x}_{i-D_{s}}^{i}$ should not be too large to capture the local contour behavior. Thus, we choose $D_{s}$ to be a fixed small number in our implementation. Some examples of $s(\mathbf{w})$ are shown in Fig. 5.

Combining the likelihood and prior terms, we get the negative log of the posterior,

[TABLE]

V Rate Term

We losslessly encode a chosen DCC string $\hat{\mathbf{x}}$ using arithmetic coding [10]. Specifically, to implement arithmetic coding using local statistics, each DCC symbol $\hat{x}_{i}\in\mathcal{A}$ is assigned a conditional probability $P(\hat{x}_{i}|\hat{\mathbf{x}}_{1}^{i-1})$ given its all previous symbols $\hat{\mathbf{x}}_{1}^{i-1}$. The rate term $R(\hat{\mathbf{x}})$ is thus approximated as the summation of negative log of conditional probabilities of all symbols in $\hat{\mathbf{x}}$:

[TABLE]

To compute $R(\hat{\mathbf{x}})$, one must assign conditional probabilities $P(\hat{x}_{i}|\hat{\mathbf{x}}_{1}^{i-1})$ for all symbols $\hat{x}_{i}$. In our implementation, we use the variable-length context tree model [20] to compute the probabilities. Specifically, to code contours in the target image, one context tree is trained using contours in a set of training images666For coding contours in the target frame of a video, the training images are the earlier coded frames. which have correlated statistics with the target image. Next we introduce the context tree model, then discuss our construction of a context tree using prediction by partial matching (PPM)777While the computation of conditional probabilities are from PPM, our construction of a context tree is novel and stands as one key contribution in this paper. [23].

V-A Definition of Context Tree

We first define notations related to the contours in the set of training images. Denote by $\mathbf{x}(m)$, $1\leq m\leq M$, the $m$-th DCC string in the training set $\mathcal{X}=\{\mathbf{x}(1),\ldots,\mathbf{x}(M)\}$, where $M$ denotes the total number of DCC strings in $\mathcal{X}$. The total number of symbols in $\mathcal{X}$ is denoted by $L=\sum^{M}_{m=1}L_{\mathbf{x}(m)}$. Denote by $\mathbf{u}\mathbf{v}$ the concatenation of sub-strings $\mathbf{u}$ and $\mathbf{v}$.

We now define $N(\mathbf{u})$ as the number of occurrences of sub-string $\mathbf{u}$ in the training set $\mathcal{X}$. $N(\mathbf{u})$ can be computed as:

[TABLE]

where $\mathbf{1}(\mathbf{c})$ is an indicator function that evaluates to $1$ if the specified binary clause $\mathbf{c}$ is true and [math] otherwise.

Denote by $P(x|\mathbf{u})$ the conditional probability of symbol $x$ occurring given its previous sub-string is $\mathbf{u}$, where $x\in\mathcal{A}$. Given training data $\mathcal{X}$, $P(x|\mathbf{u})$ can be estimated using $N(\mathbf{u})$ as done in [39],

[TABLE]

where $\delta$ is a chosen parameter for different models.

Given $\mathcal{X}$, we learn a context model to assign a conditional probability to any symbol given its previous symbols in a DCC string. Specifically, to calculate the conditional probability $P(\hat{x}_{i}|\hat{\mathbf{x}}^{i-1}_{1})$, the model determines a context $\mathbf{w}$ to calculate $P(\hat{x}_{i}|\mathbf{w})$, where $\mathbf{w}$ is a prefix of the sub-string $\hat{\mathbf{x}}^{i-1}_{1}$, i.e., $\mathbf{w}=\hat{\mathbf{x}}^{i-1}_{i-l}$ for some context length $l$:

[TABLE]

$P(\hat{x}_{i}|\mathbf{w})$ is calculated using (20) given $\mathcal{X}$. The context model determines a unique context $\mathbf{w}$ of finite length for every possible past $\hat{\mathbf{x}}^{i-1}_{1}$. The set of all mappings from $\hat{\mathbf{x}}^{i-1}_{1}$ to $\mathbf{w}$ can be represented compactly as a context tree.

Denote by $\mathcal{T}$ the context tree, where $\mathcal{T}$ is a ternary tree: each node has at most three children. The root node has an empty sub-string, and each child node has a sub-string $\mathbf{u}x$ that is a concatenation of: i) its parent’s sub-string $\mathbf{u}$ if any, and ii) the symbol $x$ (one of l, s and r) representing the link connecting the parent node and itself in $\mathcal{T}$. An example is shown in Fig. 6. The contexts of the tree $\mathcal{T}$ are the sub-strings of the context nodes—nodes that have at most two children, i.e., the end nodes and the intermediate nodes with fewer than three children. Note that $\mathcal{T}$ is completely specified by its set of context nodes and vice versa. For each $\hat{\mathbf{x}}^{i-1}_{1}$, a context $\mathbf{w}$ is obtained by traversing $\mathcal{T}$ from the root node to the deepest context node, matching symbols $\hat{x}_{i-1},\hat{x}_{i-2},\ldots$ into the past. We can then rewrite (21) as follows:

[TABLE]

V-B Construction of Context Tree by Prediction by Partial Matching (PPM)

The PPM algorithm is considered to be one of the best lossless compression algorithms [20], which is based on the context tree model. Using PPM, all the possible sub-strings $\mathbf{w}$ with non-zero occurrences in $\mathcal{X}$, i.e., $N(\mathbf{w})>0$, are contexts on the context tree. The key idea of of PPM is to deal with the zero frequency problem when estimate $P(\hat{x}_{i}|\mathbf{w})$, where sub-string $\hat{x}_{i}\mathbf{w}$ does not occur in $\mathcal{X}$, i.e., $N(\hat{x}_{i}\mathbf{w})=0$. In such case, using (20) to estimate $P(\hat{x}_{i}|\mathbf{w})$ would result in zero probability, which cannot be used for arithmetic coding. When $N(\hat{x}_{i}\mathbf{w})=0$, $P(\hat{x}_{i}|\mathbf{w})$ is estimated instead by reducing the context length by one, i.e., $P(\hat{x}_{i}|\mathbf{w}_{2}^{|\mathbf{w}|})$. If sub-string $\hat{x}_{i}\mathbf{w}_{2}^{|\mathbf{w}|}$ still does not occur in $\mathcal{X}$, the context length is further reduced until symbol $\hat{x}_{i}$ along with the shortened context occurs in $\mathcal{X}$. Let $\mathcal{A}_{\mathbf{w}}$ be an alphabet in which each symbol along with the context $\mathbf{w}$ occurs in $\mathcal{X}$, i.e., $\mathcal{A}_{\mathbf{w}}=\{x|N(x\mathbf{w})>0,x\in\mathcal{A}\}$. Based on the PPM implemented in [23], $P(\hat{x}_{i}|\mathbf{w})$ is computed using the following (recursive) equation:

[TABLE]

To construct $\mathcal{T}$, we traverse the training data $\mathcal{X}$ once to collect statistics for all potential contexts. Each node in $\mathcal{T}$, i.e., sub-string $\mathbf{u}$, has three counters which store the number of occurrences of sub-strings $l\mathbf{u}$, $s\mathbf{u}$ and $r\mathbf{u}$, i.e., $N(l\mathbf{u})$, $N(s\mathbf{u})$ and $N(r\mathbf{u})$. To reduce memory requirement, we set an upper bound $D$ on the maximum depth of $\mathcal{T}$. As done in [19], we choose the maximum depth of $\mathcal{T}$ as $D=\lceil\ln{L}/\ln{3}\rceil$, which ensures a large enough $D$ to capture natural statistics of the training data of length $L$. $\mathcal{T}$ is constructed as described in Algorithm 1.

The complexity of the algorithm is $O(D\,L)$. To estimate the code rate of symbol $\hat{x}_{i}$, we first find the matched context $\mathbf{w}$ given past $\hat{\mathbf{x}}_{i-D}^{i-1}$ by traversing the context tree $\mathcal{T}$ from the root node to the deepest node, i.e., $\mathbf{w}=\hat{\mathbf{x}}^{i-1}_{i-|\mathbf{w}|}$, and then compute the corresponding conditional probability $P(\hat{x}_{i}|\mathbf{w})$ using (23).

In summary, having defined the likelihood, prior and rate terms, our Lagrangian objective (9) can now be rewritten as:

[TABLE]

We describe a dynamic programming algorithm to minimize the objective optimally in Section VI.

VI Optimization Algorithm

Before processing the detected contours in the target image, a context tree $\mathcal{T}$ is first computed by PPM using a set of training images as discussed in Section V. The rate term $R(\hat{\mathbf{x}})$ in (24) is computed using $\mathcal{T}$. We describe a dynamic programming (DP) algorithm to minimize (24) optimally, then analyze its complexity. Finally, we design a total suffix tree (TST) to reduce the complexity.

VI-A Dynamic Programming Algorithm

To simplify the expression in the objective, we define a local cost term $f(\hat{\mathbf{x}}_{i-D}^{i})$ combining the prior and the rate term as:

[TABLE]

Note that the maximum depth $D=\lceil\ln{L}/\ln{3}\rceil$ of $\mathcal{T}$ is much larger than $D_{s}$. We assume here that the first $D$ estimated symbols $\hat{\mathbf{x}}^{D}_{1}$ are the observed $\mathbf{y}^{D}_{1}$, and that the last estimated edge is correct, i.e., $\mathbf{e}_{\hat{\mathbf{x}}}(L_{\hat{\mathbf{x}}})=\mathbf{e}_{\mathbf{y}}(L_{\mathbf{y}})$.

In summary, our algorithm works as follows. As we examine each symbol in observed $\mathbf{y}$, we identify an “optimal” state traversal through our 3-state Markov model—one that minimizes objective (24)—via two recursive functions. The optimal state traversal translates directly to an estimated DCC string $\hat{\mathbf{x}}$, which is the output of our algorithm.

Denote by $G_{i}(\hat{\mathbf{x}}_{i-D}^{i-1},\mathbf{e},j-1)$ the minimum cost for estimated $\hat{\mathbf{x}}$ from the $i$-th symbol onwards, given that we are in the good state with a set of $D$ previous symbols (context) $\hat{\mathbf{x}}^{i}_{i-D}$, and last edge is $\mathbf{e}$ which is the same as the $(j-1)$-th edge in observed $\mathbf{y}$. If we select one additional symbol $\hat{x}_{i}=y_{j}$, then we remain in the good state, incurring a local cost $f(\hat{\mathbf{x}}^{i}_{i-D})$, plus a recursive cost $G_{i+1}(\,)$ for the remaining symbols in string $\hat{\mathbf{x}}$ due to new context $\hat{\mathbf{x}}^{i}_{i-D+1}$.

If instead we choose one additional symbol $\hat{x}_{i}\neq y_{j}$, then we start a new burst error event, incurring a local cost $(c_{0}+c_{2})$ for the new event, in addition to $f(\hat{\mathbf{x}}^{i}_{i-D})$. Entering bad states, we use $B_{i+1}(\,)$ for recursive cost instead.

$B_{i}(\hat{\mathbf{x}}_{i-D}^{i-1},\mathbf{e},j-1)$ is similarly computed as $G_{i}(\hat{\mathbf{x}}_{i-D}^{i-1},\mathbf{e},j-1)$, except that if selected symbol $\hat{x}_{i}$ has no corresponding edge in $\mathbf{y}$, then we add $c_{1}$ to account for an additional error symbol. If selected symbol $\hat{x}_{i}$ has a corresponding edge in $\mathbf{y}$, then this is the end of the burst error event, and we return to good state (recursive call to $G_{i+1}(\,)$ instead). In this case, we must account for the change in length between $\hat{\mathbf{x}}$ and $\mathbf{y}$ due to this burst error event, weighted by $c_{2}$. $G_{i}(\hat{\mathbf{x}}_{i-D}^{i-1},\mathbf{e},j-1)$ and $B_{i}(\hat{\mathbf{x}}_{i-D}^{i-1},\mathbf{e},j-1)$ are defined in (26) and (27), respectively. $\mathbf{1}(\mathbf{c})$ is an indicator function as defined earlier. $\mathbf{e}=\{(m,n),d\}$ denotes the $(i-1)$-th edge of $\hat{\mathbf{x}}$, and $v(\mathbf{e},\hat{x}_{i})$ denotes the next edge given that the next symbol is $\hat{x}_{i}$. $j-1$ is the index of matched edge in $\mathbf{y}$.

Because we assume that the denoised DCC string $\hat{\mathbf{x}}$ is no longer than the observed $\mathbf{y}$, if $i>L_{\mathbf{y}}$ or $k-j<0$, we stop the recursion and return infinity to signal an invalid solution.

VI-B Complexity Analysis

The complexity of the DP algorithm is bounded by the size of the DP tables times the complexity of computing each table entry. Denote by $Q$ the number of possible edge endpoint locations. Looking at the arguments of the two recursive funcions $G_{i}(\,)$ and $B_{i}(\,)$, we see that DP table size is $O(3^{D}\,4Q\,L_{\mathbf{y}}^{2})$. To compute each entry in $B_{i}(\,)$, for each $\hat{\mathbf{x}}_{i}\in\mathcal{A}$, one must check for matching edge in $\mathbf{y}$, hence the complexity is $O(3L_{\mathbf{y}})$. Hence the total complexity of the algorithm is $O(3^{D}QL_{\mathbf{y}}^{3})$, which is polynomial time with respect to $L_{\mathbf{y}}$.

VI-C Total Suffix Tree

When the training data is large, $D$ is also large, resulting in a very large DP table size due to the exponential term $3^{D}$. In (26) and (27) when calculating local cost $f(\hat{\mathbf{x}}_{i-D}^{i})$, actually the context required to compute rate is $\mathbf{w}=\hat{\mathbf{x}}_{i-|\mathbf{w}|}^{i-1}$, where the context length $|\mathbf{w}|$ is typically smaller than $D$ because the context tree $\mathcal{T}$ of maximum depth $D$ is variable-length. Thus, if we can, at appropriate recursive calls, reduce the “history” from $\hat{\mathbf{x}}^{i}_{i+1-D}$ of length $D$ to $\hat{\mathbf{x}}^{i}_{i+1-k}$ of length $k$, $k<D$, for recursive call to $G_{i+1}()$ or $B_{i+1}$ in (26) or (27), then we can reduce the DP table size and in turn the computation complexity of the DP algorithm.

The challenge is how to retain or “remember” just enough previous symbols $\hat{x}_{i},\hat{x}_{i-1},\ldots$ during recursion so that the right context $\mathbf{w}$ can still be correctly identified to compute rate at a later recursive call. The solution to this problem can be described simply. Let $\mathbf{w}$ be a context (context node) in context tree $\mathcal{T}$. A suffix of length $k$ is the first $k$ symbols of $\mathbf{w}$, i.e., $\mathbf{w}^{k}_{1}$. Context $\mathbf{w}$ at a recursive call must be a concatenation of a chosen $i$-th symbol $\hat{x}_{i}=w_{|\mathbf{w}|}$ during a previous recursive call $G_{i}()$ or $B_{i}()$ and suffix $\mathbf{w}^{|\mathbf{w}|-1}_{1}$. It implies that suffix $\mathbf{w}^{|\mathbf{w}|-1}_{1}$ must be retained during the recursion in (26) or (27) for this concatenation to $\mathbf{w}$ to take place at a later recursion. To concatenate to suffix $\mathbf{w}^{|\mathbf{w}|-1}_{1}$ at a later recursive call, one must retain its suffix $\mathbf{w}^{|\mathbf{w}|-2}_{1}$ at an earlier call. We can thus generalize this observation and state that a necessary and sufficient condition to preserve all contexts $\mathbf{w}$ in context tree $\mathcal{T}$ is to retain all suffixes of $\mathbf{w}$ during the recursion.

All suffixes of contexts in $\mathcal{T}$ can be drawn collectively as a tree as follows. For each suffix $\mathbf{s}$, we trace $\mathbf{s}$ down from the root of a tree according to the symbols in $\mathbf{s}$, creating additional nodes if necessary. We call the resulting tree a total suffix tree (TST)888An earlier version of TST was proposed in [38] for a different contour coding application., denoted as $\mathcal{T}_{s}$. By definition, $\mathcal{T}$ is a sub-tree of $\mathcal{T}_{s}$. Further, $\mathcal{T}_{s}$ is essentially a concatenation of all sub-trees of $\mathcal{T}$ at root. Assuming $\mathcal{T}$ has $K$ contexts, each of maximum length $D$. Each context can induce $O(D)$ additional context nodes in TST $\mathcal{T}_{s}$. Hence TST $\mathcal{T}_{s}$ has $O(KD)$ context nodes.

Fig. 7 illustrates one example of TST derived from the context tree shown in Fig. 6. TST $\mathcal{T}_{s}$ can be used for compact DP table entry indexing during recursion (26) and (27) as follows. When an updated history $\hat{\mathbf{x}}_{i+1-D}^{i}$ is created from a selection of symbol $\hat{x}_{i}$, we first truncate $\hat{\mathbf{x}}_{i+1-D}^{i}$ to $\hat{\mathbf{x}}_{i+1-k}^{i}$, where $\hat{\mathbf{x}}_{i+1-k}^{i}$ is the longest matching string in $\mathcal{T}_{s}$ from root node down. The shortened history $\hat{\mathbf{x}}_{i+1-k}^{i}$ is then used as the new argument for the recursive call. Practically, it means that only DP table entries of arguments $\hat{\mathbf{x}}_{i+1-k}^{i}$ that are context nodes of TST $\mathcal{T}_{s}$ will be indexed. Considering the complexity induced by the prior term, the complexity is thus reduced from original $O(3^{D}QL_{\mathbf{y}}^{3})$ to $O((KD+3^{D_{s}})QL_{\mathbf{y}}^{3})$, which is now polynomial in $D$.

VII Experimental Results

VII-A Experimental Setup

We used two computer-generated (noiseless) depth sequences: Dude (800$\times$400) and Tsukuba (640$\times$480), and two natural (noisy) multiview silhouette sequences: Model1 (768$\times$1024) and Model2 (768$\times$1024) from Microsoft Research. 10 frames were tested for each sequence. The multiview silhouette sequences are extracted using the equipment setup in [40], where eight views are taken for one frame of each sequence. The extracted silhouettes are compressed and transmitted for further silhouette-based 3D model reconstruction [41, 40]. Fig. 8 shows one view of Model1 and Model2.

We used gradient-based edge detection [4] to detect contours. For the silhouette sequences, the detected contours were noisy at acquisition. Note that before detecting the contours for silhouette sequences, we first used a median filter999http://people.clarkson.edu/ hudsonb/courses/cs611/ to remove the noise inside and outside the silhouette. For the computer-generated depth images, we first injected noise to the depth images assuming that the pixels along the contours were corrupted by iid noise: for each edge of the contour, the corruption probability was fixed at $\delta$. If corrupted, the pixel along one side of this edge was replaced by the pixel from the other side (side was chosen with equal probability). The noisy contours were then detected from the noisy depth images. Note that for the computer-generated depth images, the ground truth contours were also detected from the original depth images. We tested two different noise probabilities, $\delta=10\%$ and $\delta=30\%$.

To code contours in a given frame of a depth / silhouette sequence, previous two frames of the same sequence were used to train the context tree. Unless specified otherwise, we set $D_{s}=4$ and $\beta=3$ in (15) for all experiments.

The quality of each decoded contour was evaluated using sum of squared distance distortion metric as proposed in [26]: the sum of squared minimum distance from the points on each decoded contour $\mathbf{\hat{x}}$ to the ground truth contour $\mathbf{x}$, denoted as $D(\hat{\mathbf{x}},\mathbf{x})$,

[TABLE]

$d(\mathbf{e}_{\hat{\mathbf{x}}}(i),\mathbf{x})$ is the minimum absolute distance between coordinate $(\hat{m},\hat{n})$ of the $i$-th edge $\mathbf{e}_{\hat{\mathbf{x}}}(i)$ of $\hat{\mathbf{x}}$ and the segment derived by all edges of $\mathbf{x}$,

[TABLE]

where $(m_{j},n_{j})$ is the coordinate of the $j$-th edge $\mathbf{e}_{\mathbf{x}}(j)$ of $\mathbf{x}$. The main motivation for considering this distortion metric stems from the popularity of the sum of squared distortion measure, which measures the aggregate errors. Note that some other distortion metrics can also be adopted according to different applications, such as the absolute area between the decoded and ground truth contours, or the depth values change due to the decoded contour shifting compared to the ground truth contour. For the silhouette sequences without ground truth contours, we set $\lambda=0$ in (24) to get the best denoising result and use it as the ground truth.

We compared performance of four different schemes. The first scheme, Gaussian-ORD, first denoised the contours using a Gaussian filter [33], then encoded the denoised contours using a vertex based lossy contour coding method [28]. The second scheme, Lossy-AEC, first denoised the noisy contours using an irregularity-detection method [4], then encoded the denoised contours using a DCC based lossy contour coding method [25]. The third scheme, Separate, first denoised the contours using our proposal by setting $\lambda=0$, then used PPM to encoded the denoised contours. The fourth scheme, Joint, is our proposal that performed joint denoising / compression of contours.

VII-B Transition Probability Estimation

For each computer-generated depth sequence, the three transition probabilities in our three-state Markov model, i.e., $p$, $q_{1}$ and $q_{2}$, were estimated from average noise statistics by comparing the noisy contours and the ground truth contours. As described in Section IV-A, we traversed the noisy DCC string $\mathbf{y}$ along with the corresponding ground truth DCC string $\mathbf{x}$ to count the number of transitions between the three states, i.e., 0, 1 and 2. The transition probability from state101010$j$ can be the same as $i$. $i$ to state $j$ is thus computed as the number of transitions from state $i$ to state $j$ divided by the number of transitions started from state $i$.

For the silhouette sequences without ground truth contours, we estimated the transition probabilities via an alternating procedure, commonly used when deriving model parameters without ground truth data in machine learning [42, 43]. Specifically, given a noisy DCC string $\mathbf{y}$, we first assigned initial values to the transition probabilities and used them to denoise $\mathbf{y}$ to get a MAP solution $\mathbf{x}_{MAP}^{\prime}$ ($\lambda=0$ in (24)). Treating $\mathbf{x}_{MAP}^{\prime}$ as ground truth and comparing $\mathbf{y}$ and $\mathbf{x}_{MAP}^{\prime}$, we updated the transition probabilities. Then, we used the updated transition probabilities to denoise $\mathbf{y}$ again to get a new MAP solution and computed the new transition probabilities by comparing $\mathbf{y}$ and the new MAP solution. This process was repeated until the transition probabilities converged.

Fig. 9 shows the changes in model parameters $p$, $q_{1}$ and $q_{2}$ across the iterations of the alternating procedure using different initial values for Model1 and Model2. We see that for each silhouette sequence with different initial values, the parameters converged to the same values in only a few iterations. Note that we can also use the alternating procedure to estimate model parameters for the computer-generated depth sequences, but the estimated results would not be as accurate as those directly computed using the available ground truth contours and the noisy contours.

VII-C Proposed Three-state Markov Model versus iid Model

To validate our proposed three-state Markov model, we compared to a naïve iid generative model for computing the likelihood. The model is as follows: for each symbol $x_{i}$ in the ground truth contour $\mathbf{x}$, we use an iid model (a coin toss) to see if it is erred in the observed contour $\mathbf{y}$. If so, we use a Poisson distribution to model the length increase in $\mathbf{y}$ over $\mathbf{x}$. Then we examine the next symbol $x_{i+1}$ and so on until the end of the DCC string. This model assumes each symbol $x_{i}$ is independent and is capable of modelling the length increase.

We adopt the widely used Akaike information criterion (AIC) [44] to compare the two models. AIC is a measure of the relative quality of statistical models for a given set of data,

[TABLE]

where $\mathcal{L}=P(\mathbf{y}|\mathbf{x})$ is the likelihood and $k$ is the number of parameters to be estimated. For the proposed three-state Markov model with three transition probabilities ($p$, $q_{1}$, $q_{2}$), $k=3$; for the i.i.d model, $k=2$ which contains $p^{\prime}$—the probability of each symbol $x_{i}$ being an error and $\lambda^{\prime}$—the average increased length per erred symbol $x_{i}$ in Poisson distribution. $p^{\prime}$ and $\lambda^{\prime}$ were estimated from average noise statistics as similar in Section VII-B. To accurately compute AIC, we only test sequence Dude and Tsukuba whose ground truth contours are available.

Table I shows the results of the negative log-likelihood and AIC of the two models. All the values of negative log-likelihood and AIC are averaged to one contour for better comparison. A smaller value of AIC represents a better fit of the model. We can see that the proposed three-state Markov model always achieved smaller values. Relative likelihood [45], computed as $exp((AIC_{prop}-AIC_{iid})/2)$, in the last column is interpreted as being proportional to the probability that the iid model is the proposed model. For example, for $Tsukuba$ with 10% noise, it means the iid model is 3.93e-2 times as probable as the proposed model to minimize the information loss. All the results show that the proposed three-state Markov model fits the data better than the iid model.

VII-D Performance of Rate Distortion

We show the rate-distortion (RD) performance of the four comparison schemes in Fig. 10 and Fig. 11. For Gaussian-ORD, RD curve was generated by adjusting the coding parameters of ORD [28]. For Lossy-AEC, RD-curve was generated by adjusting the strength of contour approximation [25]. For Joint, RD curve was obtained by varying $\lambda$. For Separate, we used the same fixed $\beta$ as in Joint to get the best denoising performance (MAP solution), then we increased $\beta$ to further smooth the contour to reduce the bit rate and obtained the RD curve.

We see that Joint achieved the best RD performance for both the depth sequences and the silhouette sequences, demonstrating the merit of our joint approach. In particular, we save about 14.45%, 36.2% and 38.33% bits on average against Separate, Lossy-AEC and Gaussian-ORD respectively for depth sequences with 10% noise, and about 15.1%, 28.27% and 32.02% for depth sequences with 30% noise, and about 18.18%, 54.17% and 28.57% respectively for silhouette sequences.

Comparing the results with 10% noise and the results with 30% noise, we find that more noise generally results in larger bit rate for Joint, Separate and Lossy-AEC. For Gaussian-ORD, the bit rate does not increase much with more noise, especially the results of Tsukuba. This was because Gaussian-ORD used a vertex-based contour coding approach; i.e., ORD encoded the locations of some selected contour pixels (vertices). After denoising using Gaussian filter, the number of selected vertices by ORD are similar for different noise probabilities, resulting in similar bit rates.

Looking at the results by Separate, when the distortion becomes larger, the bit rate reduces slowly, in some points even increases slightly, especially compared to Joint. Note that the RD curve of Separate was obtained by increasing $\beta$ to get over-smoothed denoised contour to reduce the bit rate. It shows that smoother contour dose not always ensure smaller bit rate, which in other hand verifies the importance of the rate term in the problem formulation. Thus, both the prior term and the rate term are necessary for our contour denoising / compression problem.

Compared to Lossy-AEC which encodes each DCC symbol using a limited number of candidates of conditional probabilities, our proposed Joint and Separate achieve much better coding performance by using PPM with variable length context tree model. Compared to Gaussian-ORD which encodes the contour by coding some selected vertices, our proposed Joint encodes more contour pixels, but the accurately estimated conditional probabilities enable our contour pixels to be coded much more efficiently than Gaussian-ORD.

VII-E Performance of Denoising

Fig. 12 and Fig. 13 illustrate the visual denoising results of Tsukuba and Model1. We see that the denoised contours by Joint are most visually similar to the ground truth contours. As shown inside the red circle of Fig. 12(d) (g) and Fig. 13(c), the results by Gaussian-ORD contain lots of staircase shapes. This was because that the decoded contours were constructed by connecting the selected contour pixels in order, making the result unnatural with too many staircase lines. Lossy-AEC used a pre-defined irregularity-detection approach to denoise the contours, which failed to remove the undefined noise, i.e., some noise along the diagonal directions of the contours or some noise with high irregularities, as illustrated inside the red circle of Fig. 12(e) (h) and Fig. 13(d).

VIII Conclusion

In this paper, we investigate the problem of joint denoising / compression of detected contours in images. We show theoretically that in general a joint denoising / compression approach can outperform a separate two-stage approach that first denoises then encodes the denoised contours lossily. Using a burst error model that models errors in an observed string of directional edges, we formulate a rate-constrained MAP problem to identify an optimal string for lossless encoding. The optimization is solved optimally using a DP algorithm, sped up by using a total suffix tree (TST) representation of contexts. Experimental results show that our proposal outperforms a separate scheme noticeably in RD performance.