Image Steganography using Gaussian Markov Random Field Model

TL;DR

This paper introduces a Gaussian Markov Random Field model for image steganography that improves embedding efficiency and reduces detectability by capturing pixel dependencies and optimizing KL-divergence.

Contribution

It proposes a novel GMRF-based approach with an iterative optimization scheme for secure image steganography, outperforming prior model-based methods.

Findings

GMRF model outperforms MiPOD in steganography tasks.

The approach rivals state-of-the-art HiLL in practical scenarios.

Experimental results show reduced detectability and improved embedding efficiency.

Abstract

Recent advances on adaptive steganography show that the performance of image steganographic communication can be improved by incorporating the non-additive models that capture the dependences among adjacent pixels. In this paper, a Gaussian Markov Random Field model (GMRF) with four-element cross neighborhood is proposed to characterize the interactions among local elements of cover images, and the problem of secure image steganography is formulated as the one of minimization of KL-divergence in terms of a series of low-dimensional clique structures associated with GMRF by taking advantages of the conditional independence of GMRF. The adoption of the proposed GMRF tessellates the cover image into two disjoint subimages, and an alternating iterative optimization scheme is developed to effectively embed the given payload while minimizing the total KL-divergence between cover and stego, i.e., the statistical detectability. Experimental results demonstrate that the proposed GMRF outperforms the prior arts of model based schemes, e.g., MiPOD, and rivals the state-of-the-art HiLL for practical steganography, where the selection channel knowledges are unavailable to steganalyzers.

Tables

Algorithm 1     Pseudo-code for GMRF based embedding.
Require: the ${\sigma }_{n,n}$ of a pixel, the $\rho$ associated with a two-pixel clique,
embedding payload $L$ bits, mutually disjoint cover subimages: $𝒜\cup ℬ$.
1. Initialize the change probabilities ${𝜷}^{ℬ}$ in the range $[0,0.001]$.
2. for $k=1$ to $4$ do   (GMRF with dynamical clique allocation)
3.    a) Keep the ${𝜷}^{ℬ}$ unchanged, optimize the ${𝜷}^{𝒜}$ using (32), i.e.,
${𝜷}^{𝒜}=\underset{{\beta }_{X_s}^{𝒜}}{argmin}\left\{\sum _{T_s}{D}_{KL}^{𝒜,T_s}-\lambda \left[\sum _{X_s}h\left({\beta }_{X_s}^{𝒜}\right)-L/2\right]\right\}$.
b) ${\lambda }^{𝒜}\left(k\right)=\lambda$.
4.    a) Keep the ${𝜷}^{𝒜}$ unchanged, optimize the ${𝜷}^{ℬ}$ using (34), i.e.,
${𝜷}^{ℬ}=\underset{{\beta }_{X_s}^{ℬ}}{argmin}\left\{\sum _{T_s}{D}_{KL}^{ℬ,T_s}-\lambda \left[\sum _{X_s}h\left({\beta }_{X_s}^{ℬ}\right)-L/2\right]\right\}$.
b) ${\lambda }^{ℬ}\left(k\right)=\lambda$.
5.    when $k\ge 2$
$r^{𝒜}={\lambda }^{𝒜}\left(k\right)/{\lambda }^{𝒜}\left(k-1\right)$,  $r^{ℬ}={\lambda }^{ℬ}\left(k\right)/{\lambda }^{ℬ}\left(k-1\right)$.
if ($r^{𝒜}>0.98$ && $r^{ℬ}>0.98$),    then return.
end when
6. end for
7. Compute the pixel embedding cost $d$ corresponding to ${𝜷}^{𝒜}$ and ${𝜷}^{ℬ}$:
$d=\ln \left(1/\beta -2\right).$
8. Embed $L$ bits into cover image $𝐜$ to obtain the stego image $𝐬$.

Table 2. TABLE I: Average testing error P E ¯ ¯ subscript 𝑃 𝐸 \overline{P_{E}} versus payload in bits per pixel (bpp) for MiPOD and GMRF on BOSSbase ver1.01 using ensemble 1.0 classifier with two feature sets: SRM and maxSRMd2.

Feature	Algorithm	Payload(bpp)
		0.05	0.1	0.2	0.3	0.4	0.5
SRM	MiPOD	$0.4511\pm 0.0030$	$0.4051\pm 0.0041$	$0.3274\pm 0.0026$	$0.2723\pm 0.0041$	$0.2218\pm 0.0019$	$0.1821\pm 0.0023$
	GMRF	$0.4558\pm 0.0024$	$0.4143\pm 0.0036$	$0.3428\pm 0.0038$	$0.2860\pm 0.0043$	$0.2350\pm 0.0035$	$0.1925\pm 0.0028$
maxSRMd2	MiPOD	$0.4294\pm 0.0037$	$0.3772\pm 0.0028$	$0.3037\pm 0.0028$	$0.2498\pm 0.0034$	$0.2053\pm 0.0028$	$0.1683\pm 0.0037$
	GMRF	$0.4361\pm 0.0031$	$0.3830\pm 0.0020$	$0.3078\pm 0.0036$	$0.2467\pm 0.0035$	$0.1949\pm 0.0041$	$0.1595\pm 0.0033$

Table 3. TABLE II: Average testing error P E ¯ ¯ subscript 𝑃 𝐸 \overline{P_{E}} versus payload in bits per pixel (bpp) for MiPOD with low-pass filtered costs, GMRF with low-pass filtered costs and HiLL on BOSSbase ver1.01 using ensemble 1.0 classifier with two feature sets: SRM and maxSRMd2.

Feature	Algorithm	Payload(bpp)
		0.05	0.1	0.2	0.3	0.4	0.5
SRM	HiLL	$0.4704\pm 0.0024$	$0.4330\pm 0.0029$	$0.3582\pm 0.0026$	$0.2975\pm 0.0055$	$0.2450\pm 0.0025$	$0.2018\pm 0.0027$
	MiPOD	$0.4549\pm 0.0034$	$0.4139\pm 0.0031$	$0.3417\pm 0.0031$	$0.2811\pm 0.0048$	$0.2362\pm 0.0039$	$0.1910\pm 0.0020$
	GMRF	$0.4581\pm 0.0027$	$0.4210\pm 0.0042$	$0.3530\pm 0.0031$	$0.2948\pm 0.0020$	$0.2444\pm 0.0033$	$0.2015\pm 0.0023$
maxSRMd2	HiLL	$0.4237\pm 0.0019$	$0.3732\pm 0.0045$	$0.3094\pm 0.0029$	$0.2590\pm 0.0029$	$0.2169\pm 0.0018$	$0.1789\pm 0.0033$
	MiPOD	$0.4416\pm 0.0030$	$0.3902\pm 0.0028$	$0.3231\pm 0.0036$	$0.2684\pm 0.0026$	$0.2212\pm 0.0023$	$0.1846\pm 0.0029$
	GMRF	$0.4445\pm 0.0020$	$0.3958\pm 0.0032$	$0.3240\pm 0.0020$	$0.2684\pm 0.0035$	$0.2201\pm 0.0038$	$0.1787\pm 0.0024$

Equations

$$\mathcal{N}=\{{\mathcal{N}}_k|\forall k\in \mathcal{S}\},$$

$$\begin{array}[]{l}\textbf{Positivity}:\hskip 22.76228ptp({\bf{Z}})>0,\hskip 8.5359pt\forall{\bf{Z}}\in{{\mathbb{R}}^{n}}\\ \textbf{Markovianity}:\hskip 5.69046ptp({Z_{k}}|{Z_{{{\cal S}}-\{k\}}})=p({Z_{k}}|{Z_{{{{\cal N}}_{k}}}})\end{array}$$

$$p(\mathbf{z})=((2\pi {)}^n|\mathbf{\Sigma }|{)}^{-1/2}\exp (-\frac{1}{2}(\mathbf{z}-\mathbf{\mu }{)}^T{\mathbf{\Sigma }}^{-1}(\mathbf{z}-\mathbf{\mu })),$$

$${\mathbf{X}}^m\sim N(0,{\mathbf{\Sigma }}^m),m=1,2,\cdots ,M,$$

$$\begin{array}[]{*{20}{c}}{{f}(x_{1}^{m},x_{2}^{m})=(({2\pi})^{2}{{|{{\bf{\Sigma}}^{m}}|})^{-{1\mathord{\left/{\vphantom{12}}\right.\kern-1.2pt}2}}}\exp(-\displaystyle\displaystyle\frac{1}{2}({\bf{x}}^{m})^{T}({\bf{\Sigma}}^{m})^{-1}{{\bf{x}}^{m}})},\\ {{{\bf{\Sigma}}^{m}}=\left[{\begin{array}[]{*{20}{c}}{\sigma_{1,1}^{m}}&{\sigma_{1,2}^{m}}\\ {\sigma_{2,1}^{m}}&{\sigma_{2,2}^{m}}\end{array}}\right]{\rm{}}\hskip 2.84544pt\text{and}\hskip 5.69046pt{\rm{}}({\bf{\Sigma}}^{m})^{-1}=\left[{\begin{array}[]{*{20}{c}}{\gamma_{1,1}^{m}}&{\gamma_{1,2}^{m}}\\ {\gamma_{2,1}^{m}}&{\gamma_{2,2}^{m}}\end{array}}\right]},\end{array}$$

$$F_{\Delta }(x_1^m,x_2^m)={\int }_{x_1^m-\Delta /2}^{x_1^m+\Delta /2}{\int }_{x_2^m-\Delta /2}^{x_2^m+\Delta /2}f(x_1^m,x_2^m)d{x}_2^{m}d{x}_1^m.$$

$$p_{i,j}^m={F_{\Delta }(x_1^m,x_2^m)|}_{x_1^m=i\Delta ,x_2^m=j\Delta }={\Delta }^2\cdot f(i^{\prime }\Delta ,j^{\prime }\Delta ),$$

$$\begin{array}[]{l}{p}({y_{k}}={x_{k}}+1)=\beta_{k}^{+}\\ {p}({y_{k}}={x_{k}}-1)=\beta_{k}^{-}\\ {p}({y_{k}}={x_{k}})=1-\beta_{k}^{+}-\beta_{k}^{-}\end{array},$$

$$\begin{array}[]{l}{q^{m}}({\bf{Y}}_{i,j}^{m}={\bf{X}}_{i+1,j}^{m}+{[-1,0]^{T}}|{\bf{X}}_{i+1,j}^{m})=\beta_{1}^{m}(1-2\beta_{2}^{m})\\ {q^{m}}({\bf{Y}}_{i,j}^{m}={\bf{X}}_{i,j+1}^{m}+{[0,-1]^{T}}|{\bf{X}}_{i,j+1}^{m})=(1-2\beta_{1}^{m})\beta_{2}^{m}\\ {q^{m}}({\bf{Y}}_{i,j}^{m}={\bf{X}}_{i-1,j}^{m}+{[1,0]^{T}}|{\bf{X}}_{i-1,j}^{m})=\beta_{1}^{m}(1-2\beta_{2}^{m})\\ {q^{m}}({\bf{Y}}_{i,j}^{m}={\bf{X}}_{i,j-1}^{m}+{[0,1]^{T}}|{\bf{X}}_{i,j-1}^{m})=(1-2\beta_{1}^{m})\beta_{2}^{m}\\ {q^{m}}({\bf{Y}}_{i,j}^{m}={\bf{X}}_{i-1,j-1}^{m}+{[1,1]^{T}}|{\bf{X}}_{i-1,j-1}^{m})=\beta_{1}^{m}\beta_{2}^{m}\\ {q^{m}}({\bf{Y}}_{i,j}^{m}={\bf{X}}_{i+1,j+1}^{m}+{[-1,-1]^{T}}|{\bf{X}}_{i+1,j+1}^{m})=\beta_{1}^{m}\beta_{2}^{m}\\ {q^{m}}({\bf{Y}}_{i,j}^{m}={\bf{X}}_{i+1,j-1}^{m}+{[-1,1]^{T}}|{\bf{X}}_{i+1,j-1}^{m})=\beta_{1}^{m}\beta_{2}^{m}\\ {q^{m}}({\bf{Y}}_{i,j}^{m}={\bf{X}}_{i-1,j+1}^{m}+{[1,-1]^{T}}|{\bf{X}}_{i-1,j+1}^{m})=\beta_{1}^{m}\beta_{2}^{m}\\ {q^{m}}({\bf{Y}}_{i,j}^{m}={\bf{X}}_{i,j}^{m}+{[0,0]^{T}}|{\bf{X}}_{i,j}^{m})=(1-2\beta_{1}^{m})(1-2\beta_{2}^{m})\end{array},$$

$$\begin{array}[]{l}q_{i,j}^{m}=p_{i,j}^{m}(1-2\beta_{1}^{m})(1-2\beta_{2}^{m})\\ +(p_{i-1,j}^{m}+p_{i+1,j}^{m})\beta_{1}^{m}(1-2\beta_{2}^{m})\\ +(p_{i,j-1}^{m}+p_{i,j+1}^{m})(1-2\beta_{1}^{m})\beta_{2}^{m}\\ +(p_{i-1,j-1}^{m}+p_{i-1,j+1}^{m}+p_{i+1,j-1}^{m}+p_{i+1,j+1}^{m})\beta_{1}^{m}\beta_{2}^{m}\end{array},$$

$$D(\mathbf{\beta })=D_{KL}(F_P||F_Q^{\mathbf{\beta }})\approx \frac{1}{2}{\mathbf{\beta }}^T\cdot {\nabla }^{2}D(0)\cdot \mathbf{\beta },$$

$${{\nabla }^{2}D(\mathbf{\beta })|}_{\mathbf{\beta }=0}={{\mathbf{I}}_2(\mathbf{\beta })|}_{\mathbf{\beta }=0}/\ln 2,$$

$${{\bf{I}}_{2}}({\bf{0}})=\left[{\begin{array}[]{*{20}{c}}{{I_{2}}{{({\bf{0}})}_{1,1}}}&{{I_{2}}{{({\bf{0}})}_{1,2}}}\\ {{I_{2}}{{({\bf{0}})}_{2,1}}}&{{I_{2}}{{({\bf{0}})}_{2,2}}}\end{array}}\right],$$

$$\begin{array}[]{l}{I_{2}}{({\bf{0}})_{k,l}}=E\left[{{{\left.{\left({\displaystyle\frac{{\partial\ln F_{Q}^{\bm{\beta}}({\bf{X}})}}{{\partial{\beta_{k}}}}\cdot\displaystyle\frac{{\partial\ln F_{Q}^{\bm{\beta}}({\bf{X}})}}{{\partial{\beta_{l}}}}}\right)}\right|}_{{\bm{\beta}}={\bf{0}}}}}\right],\\ \hskip 45.52458ptk\in\{1,2\},\hskip 5.69046ptl\in\{1,2\}.\end{array}$$

$$\begin{array}[]{l}{{{I}}_{2}}{({\bf{0}})_{1,1}}={\left.{\sum\limits_{i}{\sum\limits_{j}{\displaystyle\frac{{{{(\Omega_{i,j}^{1})}^{2}}}}{{{p_{i,j}}}}}}}\right|_{{\bm{\beta}}={\bf{0}}}}\\ {{{I}}_{2}}{({\bf{0}})_{1,2}}={{{I}}_{2}}{({\bf{0}})_{2,1}}={\left.{\sum\limits_{i}{\sum\limits_{j}{\displaystyle\frac{{\Omega_{i,j}^{1}\cdot\Omega_{i,j}^{2}}}{{{p_{i,j}}}}}}}\right|_{{\bm{\beta}}={\bf{0}}}}\\ {{{I}}_{2}}{({\bf{0}})_{2,2}}={\left.{\sum\limits_{i}{\sum\limits_{j}{\displaystyle\frac{{{{(\Omega_{i,j}^{2})}^{2}}}}{{{p_{i,j}}}}}}}\right|_{{\bm{\beta}}={\bf{0}}}}\end{array},$$

$$\begin{array}[]{l}\Omega_{i,j}^{1}\approx{\Delta^{4}}\cdot{\left.{\displaystyle\frac{{\partial{}^{2}{f}({x_{1}},{x_{2}})}}{{\partial{{x_{1}}^{2}}}}}\right|_{{x_{1}}=i\Delta,{x_{2}}=j\Delta}}\\ \Omega_{i,j}^{2}\approx{\Delta^{4}}\cdot{\left.{\displaystyle\frac{{\partial{}^{2}{f}({x_{1}},{x_{2}})}}{{\partial{{x_{2}}^{2}}}}}\right|_{{x_{1}}=i\Delta,{x_{2}}=j\Delta}}\end{array},$$

$$\begin{array}[]{l}\displaystyle\frac{{{\partial^{2}}{f}({x_{1}},{x_{2}})}}{{\partial{x_{1}}^{2}}}={f}({x_{1}},{x_{2}})\cdot{\left[{({\gamma_{1,1}}{x_{1}}+{\gamma_{1,2}}{x_{2}})^{2}}-{\gamma_{1,1}}\right]}\\ \displaystyle\frac{{{\partial^{2}}{f}({x_{1}},{x_{2}})}}{{\partial{x_{2}}^{2}}}={f}({x_{1}},{x_{2}})\cdot{\left[{({\gamma_{2,1}}{x_{1}}+{\gamma_{2,2}}{x_{2}})^{2}}-{\gamma_{2,2}}\right]}\end{array},$$

$$\begin{array}[]{l}{I_{2}}{{\rm{(}}{\bf{0}}{\rm{)}}_{{1,1}}}=\sum\limits_{i}{\sum\limits_{j}{\displaystyle\frac{{{{{\rm{(}}\Omega_{i,j}^{1}{\rm{)}}}^{2}}}}{{{p_{i,j}}}}}}\\ \approx{{\sum\limits_{i}{\sum\limits_{j}{\displaystyle\frac{{{{\left({{\Delta^{4}}\cdot f{\rm{(}}i\Delta,j\Delta{\rm{)}}\cdot\left[{{{{\rm{(}}{\gamma_{1,1}}i\Delta+{\gamma_{1,2}}j\Delta{\rm{)}}}^{2}}-{\gamma_{1,1}}}\right]}\right)}^{2}}}}{{{\Delta^{2}}\cdot f{\rm{(}}i\Delta,j\Delta{\rm{)}}}}}}}}\\ {\mathop{\approx}}{\Delta^{4}}\cdot\displaystyle{\iint_{{{\mathbb{R}}^{2}}}{f}({x_{1}},{x_{2}})\cdot{{\left[{{{({\gamma_{1,1}}{x_{1}}+{\gamma_{1,2}}{x_{2}})}^{2}}}-{\gamma_{1,1}}\right]}^{2}}}\mathrm{d}x_{1}\mathrm{d}x_{2},\\ \end{array}$$

$$\begin{array}[]{l}{{{I}}_{2}}{({\bf{0}})_{2,2}}=\sum\limits_{i}{\sum\limits_{j}{\displaystyle\frac{{{{(\Omega_{i,j}^{2})}^{2}}}}{{{p_{i,j}}}}}}\\ {\mathop{\approx}}{\Delta^{4}}\cdot\displaystyle{\iint_{{{\mathbb{R}}^{2}}}{f}({x_{1}},{x_{2}})\cdot{{\left[{{{({\gamma_{2,1}}{x_{1}}+{\gamma_{2,2}}{x_{2}})}^{2}}}-{\gamma_{2,2}}\right]}^{2}}}\mathrm{d}x_{1}\mathrm{d}x_{2},\\ \end{array}$$

$$\begin{array}[]{l}{{{I}}_{2}}{({\bf{0}})_{1,2}}={{{I}}_{2}}{{({\bf{0}})}_{2,1}}=\sum\limits_{i}{\sum\limits_{j}{\displaystyle\frac{{\Omega_{i,j}^{1}\cdot\Omega_{i,j}^{2}}}{{{p_{i,j}}}}}}\\ {\mathop{\approx}}{\Delta^{4}}\cdot\displaystyle{\iint_{{{\mathbb{R}}^{2}}}{f}({x_{1}},{x_{2}})\cdot{\left[{{{({\gamma_{1,1}}{x_{1}}+{\gamma_{1,2}}{x_{2}})}^{2}}}-{\gamma_{11}}\right]}}\\ \hskip 85.35826pt\cdot{{\left[{{{({\gamma_{2,1}}{x_{1}}+{\gamma_{2,2}}{x_{2}})}^{2}}}-{\gamma_{2,2}}\right]}}\mathrm{d}x_{1}\mathrm{d}x_{2}.\end{array}$$

$$\begin{array}[]{l}{\rm{}}E[X_{1}^{2}]={\sigma_{1,1}}{\rm{,}}E[{X_{1}}{X_{2}}]={\sigma_{1,2}},\\ E[X_{1}^{4}]=3{({\sigma_{1,1}})^{2}},{\rm{}}E[X_{1}^{3}{X_{2}}]=3{\sigma_{1,1}}{\sigma_{1,2}},\\ E[X_{1}^{2}X_{2}^{2}]={\sigma_{1,1}}{\sigma_{2,2}}+2{({\sigma_{1,2}})^{2}},\end{array}$$

$$\begin{array}[]{l}{I_{2}}{{\rm{(}}{\bf{0}}{\rm{)}}_{1,1}}\approx\displaystyle\frac{{2{\Delta^{4}}}}{{{{{\rm{(}}{\sigma_{1,1}}{\rm{)}}}^{2}}\cdot{{\left({1-{{{\rm{(}}{\rho_{1,{\rm{2}}}}{\rm{)}}}^{2}}}\right)}^{2}}}}\\ {I_{2}}{{\rm{(}}{\bf{0}}{\rm{)}}_{1,2}}={I_{2}}{{\rm{(}}{\bf{0}}{\rm{)}}_{2,1}}\approx\displaystyle\frac{{2{\Delta^{4}}\cdot{{{\rm{(}}{\rho_{1,{\rm{2}}}}{\rm{)}}}^{2}}}}{{{\sigma_{1,1}}\cdot{\sigma_{2,2}}\cdot{{\left({1-{{{\rm{(}}{\rho_{1,{\rm{2}}}}{\rm{)}}}^{2}}}\right)}^{2}}}}\\ {I_{2}}{{\rm{(}}{\bf{0}}{\rm{)}}_{2,2}}\approx\displaystyle\frac{{2{\Delta^{4}}}}{{{{{\rm{(}}{\sigma_{2,2}}{\rm{)}}}^{2}}\cdot{{\left({1-{{{\rm{(}}{\rho_{1,{\rm{2}}}}{\rm{)}}}^{2}}}\right)}^{2}}}}\end{array},$$

$$D_{KL}(F_P||F_Q^{\mathbf{\beta }})=\frac{(I_2{(0)}_{1,1}{{\beta }_1}^2+2I_2{(0)}_{1,2}{\beta }_1{\beta }_2+I_2{(0)}_{2,2}{{\beta }_2}^2)}{2\ln 2}.$$

$${\mathbf{r}}_n=\mathbf{G}{\mathbf{a}}_n+{\mathbf{\xi }}_n,$$

$${\hat{\mathbf{r}}}_n=\mathbf{G}({\mathbf{G}}^T\mathbf{G}{)}^{-1}{\mathbf{G}}^T{\mathbf{r}}_n.$$

$${\sigma }_{m,n}={({\mathbf{r}}_m-{\hat{\mathbf{r}}}_m)}^T({\mathbf{r}}_n-{\hat{\mathbf{r}}}_n)/(p^2-q).$$

$${\sigma }_{n,n}={({\mathbf{r}}_n-{\hat{\mathbf{r}}}_n)}^T({\mathbf{r}}_n-{\hat{\mathbf{r}}}_n)/(p^2-q).$$

$$\rho ={\sigma }_{m,n}/\sqrt{{\sigma }_{m,m}\cdot {\sigma }_{n,n}}.$$

$$\begin{array}[]{l}{{\cal A}}=\{k,|\bmod(k,2)=1\}\\ {{\cal B}}=\{k,|\bmod(k,2)=0\}\end{array},$$

$$D_{KL}^{T_s}=\frac{\sum _{c\in \mathcal{C}}([{\theta }_c]{\mathbf{\beta }}_c^T{\mathbf{I}}_{2,c}(0){\mathbf{\beta }}_c)-((\sum _{c\in \mathcal{C}}[{\theta }_c])-1)I_{1,X_s}(0){{\beta }_{X_s}}^2}{2\ln 2},$$

Full text

I Introduction

Steganography is the science and art of covert communication which aims to convey secret messages by slightly modifying the digital media (cover), e.g., images, videos and audios, to create the stego objects. To ensure the security of the covert communication, the statistical distributions of the cover source should be least possible changed for a given payload. So far, the majority of research works in steganography rely on digital images as covers, either in spatial domain for uncompressed images [1, 2, 3, 4, 5, 6, 7] or in DCT domain for JPEG compressed images [8, 9, 10, 11, 12].

At present, the mainstream of steganographic schemes for digital images is based on the concept of content adaptivity, i.e., embedding the secret data into the complex regions with rich-texture contents in images. Such adaptive steganographic schemes are generally implemented by first assigning a different embedding cost to each pixel and then embedding the secret data while minimizing the sum of costs of all changed pixels. Although the coding framework, e.g., the syndrome-trellis code (STC) [13], is well developed, which allows the steganographers to minimize an additive distortion function under the given payload constraint, most existing distortion functions are said to be heuristically defined, because one can hardly establish a direct connection between distortion and statistical detectability. Currently, several state-of-the-art content-adaptive image steganographic schemes in spatial domain are devised with such heuristically built distortion functions, which include WOW [1], S-UNIWARD [2] and HiLL [3].

MG [4] is the first steganographic scheme which provides a systematic approach for the design of distortion function based on sound mathematical principle. In specific, the distortion is associated with the steganographic Fisher Information [14], which is proportional to the Kullback-Leibler (KL) divergence between the statistical distributions of cover and stego images. The authors proposed to model the cover as a sequence of independent quantized multivariate Gaussian random variables with unequal local variances, and the embedding change probabilities for each cover pixel are derived to minimize the total KL-divergence by incorporating with the method of Lagrange multipliers for a given embedding operation (symmetric embedding) and payload. The cost values are then determined with the obtained change probabilities. Although a rather simple variance estimator was employed, the authors showed that the security performance of MG was comparable to the previous state-of-the-art HUGO algorithm [5].

The authors extended their work of MG later in [6] by replacing the multivariate Gaussian model with the generalized multivariate Gaussian (hence its name MVGG) and incorporating an improved variance estimator. The MVGG was known as model driven in the sense that it allows pentary embedding to hide large payload in complex and texture regions of cover images with a thicker-tail MVGG model. The MVGG showed comparable performance with pentary coded S-UNIWARD and HiLL on SRM [15] and maxSRMd2 [16] feature sets.

More recently, following the framework of model-based steganography, an alternative approach [7] for the design of distortion function was proposed by minimizing the power of optimal detector (MiPOD). Although MiPOD was a continuation of MVGG, it still shed some new insights into steganography design. Firstly, a closed-form expression for the power of the most powerful detector of adaptive embedding with LSB matching is derived based on a multivariate Gaussian cover image model developed in [4]. Secondly, the closed-form expression for statistical detectability allows ones to design the so-called “detectability-limited sender” that controls the size of secure payload for a given image to not exceed a target detectability level. Finally, the MiPOD derives the embedding change probabilities (selection channels) for steganography by minimizing the power of optimal detector rather than the KL divergence of the statistical distributions between cover and stego objects as employed in MG and MVGG. Equipped with the improved variance estimator for model parameters, the MiPOD can rival the state-of-the-art steganographic schemes in spatial domain, e.g., S-UNIWARD and HiLL.

In this paper, we go a step further to investigate model that can capture dependencies among spatially adjacent pixels for efficient image steganography with symmetric embedding in spatial domain. In specific, a Gaussian Markov Random Field model (GMRF) with four-element cross neighborhood is proposed to characterize the cover image, which tessellates the cover image set ${\bf{c}}$ into two disjoint subsets (sublattices) ${\bf{c}}_{e}$ and ${\bf{c}}_{o}$. For the given GMRF model, each pixel is associated with four two-pixel cliques (two horizontally and two vertically neighboring pixels), which makes up a subtree with five pixels. And the involved adjacent pixel pair is modeled as the jointly Gaussian random variables. To embed $L$ bits message in cover $\bf{c}$, we assign $L/2$ bits into sublattices ${\bf{c}}_{e}$ and ${\bf{c}}_{o}$ to generate stegos ${\bf{s}}_{e}$ and ${\bf{s}}_{o}$, respectively. The Markov Random Field (MRF) modeling is appealing because the subtrees associated with one sublattice, say ${\bf{c}}_{e}$, are conditionally independent when ${\bf{c}}_{o}$ is determined, which allows us to formulate the total KL-divergence between ${\bf{c}}$ and ${\bf{s}}$ as the sum of KL-divergence of all subtrees associated with ${\bf{c}}_{e}$ and ${\bf{c}}_{o}$ in terms of the ones of 4 neighboring pixel pairs. An alternating iterative strategy among two sublattices is applied to derive the embedding change probabilities (selection channels) and then the distortion costs for each pixel in one sublattice when the embedding payload and the selection channels for another sublattice are given by minimizing the total KL-divergence between cover and stego in the corresponding sublattice. Extensive experiments are carried out to verify the effectiveness of the proposed method (known as GMRF) using steganalyzers with rich model, e.g., SRM and maxSRMd2, on BOSSbase database [17]. The proposed GMRF outperforms its prior art with independent Gaussian model, i.e., MiPOD, in terms of secure embedding capacity for SRM, and can rival the current state-of-the-art one, i.e., HiLL, with symmetric embedding.

The rest of this paper is structured as follows. In Section II, we first introduce the notations and conventions adopted in this paper and review the necessary preliminaries on MRF. Then the cover and stego image models in terms of a pair of adjacent pixels (clique for GMRF with four-element cross neighborhood) based on jointly Gaussian distribution are derived. The GMRF model is then employed in Section III to formulate the problem of image steganography as the one of minimizing the KL-divergence between cover and stego on two disjoint sublattices, which are followed by the experimental results and analysis in Section IV. And finally the conclusion remarks are drawn in Section V.

II Notations and Image Modelling

II-A Notations and conventions

Let us first give the notations and conventions used throughout the paper. The calligraphic fonts will be used solely for sets (occasionally, other symbols may also be used for sets in the interest of ease presentation), random variables will be typeset in capital letters, while their corresponding realizations will be in lower cases. Vectors and matrices will always typeset in bold fonts. Cover and stego images ${\bf{c}}=({c_{k}})$ and ${\bf{s}}=({s_{k}})$ are addressed by their one-dimensional index set ${{\cal S}}=\{1,2,\cdots,K\}$, where ${c_{k}}$ and ${s_{k}}$ $\in\left\{{0,1,\cdots,G-1}\right\}$ and $G=256$ for 8-bit grayscale images. We also use the Iverson bracket, $\left[P\right]$, defined as $\left[P\right]$=1 when the statement $P$ is true and zero otherwise.

II-B Gaussian Markov Random Field (GMRF)

Let ${\cal S}$ be a set of $n$ discrete sites, i.e., ${{\cal S}}=\{1,2,\cdots,n\}$, and ${\bf{Z}}=({Z_{1}},{Z_{2}},\cdots,{Z_{n}})$ be a family of random variables defined on set ${\cal S}$, the neighborhood system for $\bf{Z}$ is defined as

[TABLE]

where ${\cal N}_{k}$ is the collection of sites neighboring to $k$ for which $k\notin{\cal N}_{k}$ and $k\in{{{\cal N}}_{l}}\Leftrightarrow l\in{{{\cal N}}_{k}}$. In MRF, the ordering of sites are not important, and their relationship is determined by the corresponding neighborhood system ${\cal N}$. A subset ${c}\in{{\cal S}}$ is known as a clique if each pair of different elements from $c$ is neighbor. In this paper, unless stated otherwise, we only consider the four-element cross neighborhood and the associated two-pixel cliques of horizontally or vertically neighboring pixels as shown in Fig. 1.

For the family of random variables $\bf{Z}$, each random variable $Z_{k}$ takes the value ${z_{k}}\in{\mathbb{R}}$ with probability $p({Z_{k}}={z_{k}})$, and the joint event $({Z_{1}}={z_{1}},\cdots,{Z_{n}}={z_{n}})$ is abbreviated as $\bf{Z}=\bf{z}$ with the probability $p(\bf{Z}=\bf{z})$. $\bf{Z}$ is known as a Markov Random Field (MRF) on ${\cal S}$ with respect to the neighborhood system ${\cal N}$ if and only if (iff) the following two conditions are satisfied [18]:

[TABLE]

The Markovianity says that the probability of a local event at $k$ conditioned on all the remaining events is equivalent to that conditioned on the events at the neighbors of $k$. A random vector ${\bf{Z}}={({Z_{1}},{Z_{2}},\cdots,{Z_{n}})^{T}}\in{{\mathbb{R}}^{n}}$ is called a Gaussian MRF (GMRF) w.r.t. a labeled graph ${{\cal G}}=({{\cal V}},{{\cal E}})$ (${\cal V}$: Vertices; ${\cal E}$: Edges) with mean ${\bm{\mu}}$ and precision matrix ${\bm{\Sigma}}^{-1}>0$ , iff it is Gaussian distributed with the form [19]:

[TABLE]

and ${\bm{\Sigma}}_{k,l}^{-1}\neq 0\Leftrightarrow\left\{{k,l}\right\}\in{\rm{{\cal E}}}$ for all $k\neq l$.

Once the condition for Positivity is satisfied, which is always true for practical applications, the joint probability $p({\bf{z}})$ is uniquely determined by its local conditional probabilities [20]. In other words, the $p({\bf{z}})$ is characterized with the distributions of all involved cliques associated with the given neighborhood system. We then proceed to develop the cover and stego image model based on the distributions of horizontally or vertically neighboring two-pixel clique in the next two subsections.

II-C Cover image model

As will be discussed later in Section III, a GMRF model with four-element cross neighborhood system is utilized to model the cover images in this paper. For a given uniform scalar quantizer ${Q_{\Delta}}$ with quantization step ${\Delta}$, the cover is quantized to a zero mean jointly distributed Gaussian ${Q_{\Delta}}(N({\bf{0}},{\bf{\Sigma}}))$. Without loss of generality, in practice, the mean of the image cover is removed from each pixel, and the ${\Delta}$ is generally set to be 1, in the interest of simplicity. According to [21], for a multivariate Gaussian random variable, its marginal variables are also Gaussian distributed, as a result, both the single pixel and the horizontally or vertically neighboring pixel pair in the cover images, which are associated with the adopted neighborhood system, are assumed to be Gaussian and jointly Gaussian distributed, respectively. Let ${\bf{X}}=\{{{\bf{X}}^{1}},{{\bf{X}}^{2}},\cdots,{{\bf{X}}^{M}}\}$ be the sequence of all involved two-pixel cliques in the cover image, for the $m^{th}$ clique ${{\bf{X}}^{m}}={\left[{X_{1}^{m},X_{2}^{m}}\right]^{T}}$, we have

[TABLE]

that is

[TABLE]

where ${{\bf{\Sigma}}^{m}}$ and $({\bf{\Sigma}}^{m})^{-1}$ are the covariance and precision matrices, resepctively. The corresponding Probability Mass Function (p.m.f) for the $m^{th}$ clique can be computed under the fine quantization step $\Delta$,

[TABLE]

Let $x_{1}^{m}=i\Delta,x_{2}^{m}=j\Delta$, the p.m.f $p_{i,j}^{m}$ for the $m^{th}$ clique can then be evaluated by incorporating the Mean Value Theorem (MVT):

[TABLE]

where $i^{\prime}\in(i-0.5,i+0.5)$ and $j^{\prime}\in(j-0.5,j+0.5)$.

II-D Stego image model

Unlike the MG [4] and MiPOD [7] with mutually inde-pendent embedding, we take into account the interactions among embedding changes of neighboring pixels when modifying the cover ${\bf{x}}=({x_{k}})$ to stego ${\bf{y}}=({y_{k}})$, where ${x_{k}}$ and ${y_{k}}$ are the $k^{th}$ cover and stego pixels, respectively. However, to simplify the problem, the same symmetric ternary embedding model as shown in Fig. 2 is adopted throughout the paper, i.e.,

[TABLE]

where the change probability $\beta_{k}^{+}=\beta_{k}^{-}={\beta_{k}}$, ${\beta_{k}}\in[0,{1\mathord{\left/{\vphantom{13}}\right.\kern-1.2pt}3}]$ is the change probability for the $k^{th}$ pixel of cover $\bf{x}$.

With the symmetric embedding model, the two-pixel cover clique sequence ${\bf{X}}=\{{{\bf{X}}^{m}}\}$ with p.m.f ${p^{m}}(x_{1}^{m},x_{2}^{m})$ is changed to stego clique sequence ${\bf{Y}}=\{{{\bf{Y}}^{m}}\}$ with p.m.f ${q^{m}}(y_{1}^{m},y_{2}^{m})$, where $m=1,2\cdots,M$. Let $\beta_{1}^{m}$ and $\beta_{2}^{m}$ be the change probabilities for pixels $X_{1}^{m}$ and $X_{2}^{m}$ in the $m^{th}$ clique, and ${\bf{X}}_{i,j}^{m}$ and ${\bf{Y}}_{i,j}^{m}$ denote the $m^{th}$ two-pixel cover and stego cliques which take the value $[i,j]$, respectively. Under ternary embedding with symmetric changes, there are totally 9 embedding change types as shown in Fig. 3,

[TABLE]

where $\beta_{1}^{m}$ and $\beta_{2}^{m}$ $\in\left[{0,{1\mathord{\left/{\vphantom{13}}\right.\kern-1.2pt}3}}\right]$. Therefore, the p.m.f for the $m^{th}$ clique in stego image can be written as

[TABLE]

where $p_{i,j}^{m}$ and $q_{i,j}^{m}$ are the probability values for the $m^{th}$ two-pixel cover and stego cliques which take the value $[i,j]$, respectively.

III Optimize the Statistical Undetectability by Minimizing the KL-Divergence

In this Section, we formulate the problem of secure image steganography as the one of minimization of KL-divergence between cover and stego images in terms of a series of low-dimensional clique structures by utilizing the conditional independence of GRMF. The KL-divergence for a two-pixel clique is firstly derived. And the parameter estimation associated with the GMRF is then discussed. Based on these, the total KL-divergence model between cover and stego images is developed and an alternating iterative optimization scheme is proposed to solve the model for the minimization of KL-divergence, i.e., least statistical detectability.

III-A The KL-Divergence of a two-pixel clique

For a two-pixel cover clique ${\bf{X}}={[{X_{1}},{X_{2}}]^{T}}$, the embedding modifies the two pixels with probabilities ${\bm{\beta}}={\left[{{\beta_{1}},{\beta_{2}}}\right]^{T}}$, changing to the stego clique ${\bf{Y}}={[{Y_{1}},{Y_{2}}]^{T}}$. Note that we omit the superscript $``m"$ for simplicity. For small change probabilities ${\bm{\beta}}$, the KL-divergence between cover and stego cliques can be well approximated with its leading quadratic terms (detailed derivations see Appendix A-1):

[TABLE]

where ${F_{P}}$ and $F_{Q}^{\bm{\beta}}$ are the p.m.f of the two-pixel cover clique and stego clique, respectively, ${\nabla^{2}}D({\bf{0}})$ is the second order partial derivative of $D$ w.r.t ${\bm{\beta}}$ when ${\bm{\beta}}=\bf{0}$. Similar to the result in [22], ${\nabla^{2}}D({\bm{\beta}})$ is also proportional to the Fisher information matrix ${{\bf{I}}_{2}}({\bm{\beta}})$ at ${\bm{\beta}}={\bf 0}$ (for proof see Appendix A-1), i.e.,

[TABLE]

where ${{\bf{I}}_{2}}({\bf{0}})$ is the binary steganographic Fisher Information Matrix (FIM)

[TABLE]

where

[TABLE]

Furthermore, by incorporating (10) and (7) for the derivation of (14) through Taylor expansion, we can derive the closed-form expression for FIM ${{\bf{I}}_{2}}({\bf{0}})$ (refer to Appendix A-2 for detailed derivations):

[TABLE]

where $\Omega_{i,j}^{1}=({p_{i-1,j}}+{p_{i+1,j}}-2{p_{i,j}})$, $\Omega_{i,j}^{2}=({p_{i,j-1}}+{p_{i,j+1}}-2{p_{i,j}})$, and by applying the result in (7), for fine quantization step $\Delta$, $f(i^{\prime}\Delta,j^{\prime}\Delta)\approx f(i\Delta,j\Delta)$, where $i^{\prime}\in(i-0.5,i+0.5)$, $j^{\prime}\in(j-0.5,j+0.5)$, thus we further have

[TABLE]

where $f(x_{1},x_{2})$ is the bivariate Gaussian p.d.f as shown in (5). It is observed that $\Omega_{i,j}^{1}$ and $\Omega_{i,j}^{2}$ are linearly proportional to the second order partial derivatives of ${f}({x_{1}},{x_{2}})$ (bivariate Gaussian for two-pixel clique) w.r.t. ${x_{1}}$ and ${x_{2}}$ at ${x_{1}}=i\Delta$ and ${x_{2}}=j\Delta$, respectively. On the other hand, the second order partial derivatives of ${f}({x_{1}},{x_{2}})$ in (16) can be easily obtained according to (5), i.e.,

[TABLE]

where ${\gamma_{1,1}}$, ${\gamma_{1,2}}$, ${\gamma_{2,1}}$ and ${\gamma_{2,2}}$ are the components of the precision matrix ${{\bm{\Sigma}}^{-1}}$ in (5). Based on these, substituting (17) into (16) at ${x_{1}}=i\Delta,{x_{2}}=j\Delta$, then we can have the FIM ${\bf{I}}_{2}(\bf{0})$ (${\lim\Delta\to 0}$)

[TABLE]

[TABLE]

[TABLE]

It is readily seen that the elements of FIM in (13) can be represented in terms of ${2^{th}}$ order and ${4^{th}}$ order moments of bivariate Gaussian $\bf{X}$, which, according to Isserlis’ theorem [23], can be represented as

[TABLE]

where ${\sigma_{1,2}}$ is the covariance of $X_{1}$ and $X_{2}$, ${\sigma_{1,1}}$ and ${\sigma_{2,2}}$ are the variance of $X_{1}$ and $X_{2}$, respectively. Therefore, we finally obtain the FIM ${{\bf{I}}_{2}}({\bf 0})$ as follows:

[TABLE]

where ${\rho_{1,2}}={{{\sigma_{1,2}}}\mathord{\left/{\vphantom{{{\sigma_{1,2}}}{\sqrt{{\sigma_{1,1}}{\sigma_{2,2}}}}}}\right.\kern-1.2pt}{\sqrt{{\sigma_{1,1}}\cdot{\sigma_{2,2}}}}}$ is the correlation coefficient of $X_{1}$ and $X_{2}$.

Thus far, we can rewrite the KL-divergence between a pair of two-pixel cover and stego cliques by taking advantage of (11), (12), (13) and (22),

[TABLE]

It is observed that the KL-divergence of a two-pixel clique is relevant to the correlation coefficient, variances and change probabilities of pixels in the clique, among which, the correlation coefficients for each neighboring pixel pair in the involved cliques are the predominant parameters to be determined for the incorporated Gaussian MRF model. In general, the performance of the proposed MRF based steganographic scheme is heavily dependent on the estimation accuracy of the model parameters. To this end, we follow in spirit the scheme in MiPOD [7] to estimate the pixel variance and the correlation coefficient of neighboring pixels in a clique, which will be illustrated later in next subsection.

III-B The estimation of pixel variance and correlation coefficient of neighboring pixels

In [7], Sedighi et al. proposed an elegant pixel variance estimator to build the underlying multivariate independent Gaussian model, which leads to notable performance improvement over its predecessor MG [4]. The variance estimator consists of two steps, that is: 1) suppress the image content using a denoising filter to obtain the residual image; 2) fit a local parametric model to the neighbors of each residual to obtain its variance estimation.

In light of its effectiveness, the same strategy is also generalized to develop the covariance estimator for the involved Gaussian MRF model. For the 8-bit grayscale cover image with original pixel values ${\bf{c}}=({c_{1}},{c_{2}},\cdots,{c_{n}})$, we obtain the residual ${\bf{r}}$ of the cover ${\bf{c}}$ using a two-dimensional wiener filter $F$: ${\bf{r}}={\bf{c}}-F({\bf{c}})$. Then we utilize a local parametric model to estimate the pixel variance by blockwise Maximum Likelihood Estimation (MLE) [24]. For the ${n^{th}}$ residual ${{{r}}_{n}}$ in cover ${\bf{c}}$, we fit the local parametric model [25] to the $p\times p$ neighbors of ${{{r}}_{n}}$ and model the residual expectation within the associated $p\times p$ block as follows:

[TABLE]

where ${{\bf{r}}_{n}}$ denotes the residual values inside the $p\times p$ block surrounding the ${n^{th}}$ residual ${{{r}}_{n}}$, which is a column vector of size ${p^{2}}\times 1$, ${\bf{G}}$ is the matrix of size ${p^{2}}\times q$ which defines the adopted parametric model, ${{\bf{a}}_{n}}$ is the parameter vector of size $q\times 1$, and ${{\bm{\xi}}_{n}}$ is the noise vector of the model with the variance and covariance to be estimated. Under the assumption of Gaussian noise, the estimation of the residual ${{\bf{r}}_{n}}$ is

[TABLE]

In general, we could assume the residuals within the $p\times p$ block centered on ${r_{n}}$ exhibit the similar statistical characteristics, i.e., they could be regarded as the multiple samples of ${r_{n}}$. Therefore, the covariance of neighboring pixels ${{{c}}_{m}}$ and ${{{c}}_{n}}$ in a clique can be well estimated as

[TABLE]

Similarly, we have the variance estimation for pixel ${{{c}}_{n}}$, i.e.,

[TABLE]

Given the variance and covariance estimations for pixels ${{{c}}_{m}}$ and ${{{c}}_{n}}$, we can readily obtain the correlation coefficient between them, i.e.,

[TABLE]

Considering the numerical stability and computational efficiency, we set adequate bounds for $\rho$ and ${\sigma_{n,n}}$, i.e., $\rho=min(\rho,0.99*sign(\rho))$ and ${\sigma_{n,n}}=max(0.01,{\sigma_{n,n}})$.

III-C Minimize the total KL-divergence between cover and stego by incorporating the Gaussian MRF

With the formulation of KL-divergence of a two-pixel clique in (23), we then further derive the total KL-divergence between cover and stego by incorporating the given Gaussian MRF with four-element cross neighborhood as shown in Fig. 1, where the image of size $M\times N$ is decomposed into two interleaved subimages (sublattices) ${\cal A}$ (“$+$”) and ${\cal B}$ (“$\times$”), i.e.,

[TABLE]

where $1\leq k\leq M\times N$.

For any pixel ${X_{s}}$ in subimage ${\cal A}$ (the same for ${\cal B}$), there are four cliques associated with it, i.e., ${{\bf{X}}_{{s},{t}}}=[{X_{s}},{X_{t}}]^{T}$, $t=1,\cdots,4$, which constitute a 4-ary clique tree denoted by ${T_{s}}=\{X_{s},X_{1},X_{2},X_{3},X_{4}\}$, as shown in Fig. 1. With given ${X_{s}}$ (or ${\beta_{s}}$) in ${\cal A}$, the other four pixels (${X_{t}}$) in cliques ${{\bf{X}}_{{s},{t}}}$ ($t=1,\cdots,4$) of ${T_{s}}$ are conditionally independent under the proposed Gaussian MRF model with four-element cross neighborhood. We then obtain the KL-divergence for clique tree ${T_{s}}$ (refer to Appendix B for details)

[TABLE]

where ${{\cal C}}=\left\{{{\rm{\{}}{X_{s}},{X_{t}}{\rm{\}}},|t=1,\cdots 4}\right\}$, ${{\bm{\beta}}_{c}}={[{\beta_{X_{s}}},{\beta_{X_{t}}}]^{T}}$, ${{\bf{I}}_{2,{c}}}({\bf{0}})$ is binary FIM for clique ${c}$, $I_{1,{{X_{s}}}}(0)$ is the Fisher Information for pixel $X_{s}$, and $[{\theta_{c}}]=1$ when clique ${c}$ is included in clique tree $T_{s}$ and zero otherwise.

Note that, with the underlying GMRF, for given change probabilities for pixels in ${\cal B}$, the pixels in ${\cal A}$ and their associated 4-ary clique trees are mutually independent. And so are the pixels in ${\cal B}$. Let the total payload be $L$ bits, we assign half of the payload to ${\cal A}$ and ${\cal B}$, respectively. Without loss of generality, we take the embedding in ${\cal A}$ for example, the allowable payload that can be embedded into subimage ${\cal A}$ is the sum of entropies of $\left\{{\beta_{X_{s}}^{{\cal A}},\beta_{X_{s}}^{{\cal A}},1-2\beta_{X_{s}}^{{\cal A}}}\right\}$, i.e.,

[TABLE]

where $\beta_{X_{s}}^{{\cal A}}$ is the change probability for the ${s^{th}}$ pixel in ${\cal A}$, and $h(x)=-2x\log x-(1-2x)\log(1-2x)$ is expressed in bits. The minimum distortion embedding for ${\cal A}$ in terms of $\beta_{X_{s}}^{{\cal A}}$ is then formulated as the minimization of total KL-divergence between cover and stego of subimage ${\cal A}$ subject to the payload constraint (31), which can be solved using the method of Lagrange multipliers,

[TABLE]

where $D_{KL}^{{{\cal A}},T_{s}}$ is the KL-divergence expressed in (30) for the ${s^{th}}$ clique tree in ${\cal A}$. Differentiating the objective w.r.t. $\beta_{X_{s}}^{{\cal A}}$ gives (refer to Appendix C):

[TABLE]

where ${\Gamma_{X_{s}}}=\sum\limits_{{c}\in{{\cal C}}}{[\theta_{c}^{T_{s}}]{\bf{I}}_{2,{c}}^{T_{s}}{{({\bf{0}})}_{2,2}}-(}(\sum\limits_{{c}\in{{\cal C}}}{[\theta_{c}^{T_{s}}]})-1)I_{1,{X_{s}}}^{T_{s}}(0)$ and ${\Lambda_{X_{s}}}=\sum\limits_{{c}\in{{\cal C}}}{[\theta_{c}^{T_{s}}]{\bf{I}}_{2,{c}}^{T_{s}}{{({\bf{0}})}_{1,2}}\beta_{X_{s\_c}}^{{\cal B}}}$, $\beta_{X_{s\_c}}^{{\cal B}}$ is the change probability for pixel ${X_{s\_c}}$ corresponding to the clique $c=\{{X_{s}},{X_{s\_c}}\}{\rm{}}$ of the 4-ary clique tree $T_{s}$, which is located in subimage $\cal B$ and is fixed in the optimization process, ${\bf{I}}_{2,{c}}^{T_{s}}({\bf{0}})$ and $I_{1,{X_{s}}}^{T_{s}}(0)$ are the binary FIM for clique ${c}$ and FI associated with the ${s^{th}}$ pixel in ${\cal A}$, respectively. The algorithm in [26] is then adopted to obtain $\bm{\beta}^{{\cal A}}$.

Similarly, with the solved $\bm{\beta}^{{\cal A}}$, we have the optimal embedding for subimage ${\cal B}$ to obtain $\bm{\beta}^{{\cal B}}$

[TABLE]

Therefore, an alternating iterative optimization scheme is developed to embed $L$ bits in cover image while minimizing the total KL-divergence by incorporating the proposed GMRF model as shown in Fig. 4. In specific, randomly initialize the $\bm{\beta}^{{\cal B}}$ in subimage ${\cal B}$ and keep them unchanged, solve (32) for $\bm{\beta}^{{\cal A}}$ in ${\cal A}$. Update the change probabilities for pixels in ${\cal A}$ with the newly obtained $\bm{\beta}^{{\cal A}}$ and then keep them unchanged, solve (34) for $\bm{\beta}^{{\cal B}}$ in ${\cal B}$. The process continues until it converges to obtain the optimal $\bm{\beta}^{{\cal A}}$ and $\bm{\beta}^{{\cal B}}$ for embedding payload of $L$ bits into cover image with minimum distortion. It is noted that: 1) the update of $\bm{\beta}^{{\cal A}}$ is affected somehow by the neighboring $\bm{\beta}^{{\cal B}}$ through the underlying GMRF model; 2) the problems specified in (32) and (34) are the ones of convex optimization, which are bound to converge to the globally optimal solutions. The readers are advised to refer to the pseudo-code (Algorithm 1) with practical consideration to better understand the process.

We then proceed to consider the practical issues in implementing the GMRF based image steganography. Recall that the interactions between neighboring pixels are explored by taking advantages of the proposed GMRF with four-element cross neighborhood, which doesn’t go well enough for embedding in rich texture regions of cover images. Fig. 5 shows a portion of image (1013.pgm) from BOSSbase [17], where the red boxes correspond to the texture regions and the black and white points inside images represent the ’$-$1’ and ’$+$1’ embeddings, respectively. The paradigm of adaptive image steganography encourages more data to be embedded in the texture regions of images by assigning low costs to changes to regions which exhibit less correlations among neighboring pixels. The red boxes in Fig. 5 (c) shows the embedding distribution using MiPOD [7] with a multivariate independent Gaussian model, which is appropriate to characterize the rich texture regions of images, while Fig. 5 (d) shows slightly sparser distributions in the same regions with the proposed GMRF, where all the four two-pixel cliques are included in the clique tree. This is because, for the two less correlated pixels in a clique, say, $X_{1}$ (low cost) and $X_{2}$ (high cost), the embedding cost for $X_{1}$ would increase during the alternating iterative optimization process, which may lead to less embeddings in the texture regions. Therefore, if two horizontally/vertically neighboring pixels are less correlated to some extents, the relevant clique should be disconnected from the GMRF model. Note that the correlation between two neighboring pixels is in connection with their variances, the involved two pixels are less correlated as long as one or two of them have relatively large variances. On the other hand, the ultimate goal for GMRF based optimization is to find the optimal change probability $\beta$ of each pixel $x$ in the cover image for minimum distortion embedding, while $\beta$ is directly proportional to pixel’s variance, thus the change probability $\beta$ can be used as an effective measure to dynamically allocate the relevant cliques in the alternating iterative optimization process, i.e., let ${{\bf{X}}_{c}}={[{X_{s}},{X_{t}}]^{T}}$ be the clique ${c}=\{{X_{s}},{X_{t}}\}$ ($t=1,..,4$) associated with the ${s^{th}}$ clique tree, we have

[TABLE]

where $\theta_{c}^{T_{s}}$ is the statement to determine if clique ${{\bf{X}}_{c}}$ is included in the clique tree $T_{s}$, ${\beta_{T}}$ is the threshold, ${\beta_{X_{s}}}$ and ${\beta_{X_{t}}}$ are the change probabilities for the two pixels in clique ${{\bf{X}}_{c}}$, respectively. The introduction of (35) can then be regarded as the dynamical configuration of the initial neighborhood system ${\cal N}$ for GMRF based embedding with practical consideration. Fig. 5 (e) and (f) show the embedding distributions of GMRF in texture regions when the dynamical allocation scheme is applied and all the cliques are excluded in the clique trees, respectively. It is observed that the GMRF has denser embedding densities than MiPOD in red blocks, while shows similar densities with MiPOD in highly-textured regions. Note that the red blocks in Fig. 5 may consists of highly and medium textured regions, and the performance gains of the proposed GMRF over independent Gaussian model based schemes are mainly from the medium and less textured regions in images.

As a summary, we finally give the pseudo-code (Algorithm 1) for GMRF based embedding with practical consideration.

IV Experimental Results and Analysis

IV-A Experiment setups

All the experiments in this Section are carried out on image database BOSSbase ver1.01 [17] which contains 10,000 gray-scale images of size $512\times 512\times 8$ bits. All tested schemes are simulated at their corresponding payload-distortion bound, which execute the embedding modifications with the probabilities. The block size for variances and correlation coefficients estimation is $9\times 9$, and the threshold for dynamical clique allocation is set as ${\beta_{T}}=0.1$. Two state-of-the-art feature sets for spatial images, i.e., SRM [15] and its selection-channel-aware version maxSRMd2 [16] are used to evaluate the empirical security performance of the proposed GMRF based scheme and other competing methods. The Fisher Linear Discriminant ensemble [27] is also adopted in our experiments to train the binary classifier. Half of the cover and stego images will be used as the training set for the ensemble classifiers, and the remaining half will be used as test set to evaluate the trained classifier. And the security performance is quantified as the minimal total probability of error under equal priors achieved on the test set by ten times of randomly testing, denoted as ${\overline{P_{E}}}$.

IV-B Performance comparison with other competing scheme

We compare the proposed GMRF with the baseline of the state-of-the-art model based scheme MiPOD without smoothing on its Fisher information. Table I show the security performance of the two schemes against SRM and maxSRMd2. It is observed that the GMRF consistently outperforms MiPOD by a clear margin for SRM across the tested payloads, indicating that the proposed GMRF can better preserve the statistical distribution of stego images after embedding and lead to less detectibility. For the selection-channel-aware maxSRMd2, however, our GMRF has slightly better performance than MiPOD for low payloads ($R\leq 0.2{\rm{bpp}}$), and tends to be inferior to MiPOD for medium and large payloads ($R\geq 0.3{\rm{bpp}}$), which is the range of insecurity (the corresponding $\overline{P_{E}}<0.4$). This is most likely because the GMRF could better characterize the images than MiPOD, which can be acquired by maxSRMd2 through the change probabilities $\beta$ to more effectively detect the stego images.

IV-C The effect of smoothing operations on GMRF and MiPOD

Nowadays, it is a common practice to boost the security performance of steganographic schemes by smoothing the embedding costs using a low-pass filter [28]. This can be explained from the perspective of maximum entropy principle for discrete source in information theory that the smoothing operation on embedding costs tends to uniformize the change probabilities in local regions and thus increase the embedding entropy in highly textured regions. In addition, the smoothing operation can also spread the high costs of pixels into their neighborhood which would make the embedding more conservative in textured edges. In [7], the authors take a similar measure to improve the performance of MiPOD by smoothing the Fisher information, which is closely correlated with the embedding cost. While in the paper, for ensuring the fairness of the comparative experiment, we will choose to smooth the embedding cost of MiPOD directly, since smoothing the binary Fisher information in FIM for GMRF is inappropriate.

In our implementation, we adopt a low-pass filter with support as the one used in MiPOD, i.e., a $7\times 7$ average filter. Table II show the performance comparison of HiLL, and low-pass filtered GMRF and MiPOD. It is observed that, for SRM, the proposed GMRF with low-pass filtered cost outperforms GMRF without low-pass filtering and MiPOD with low-pass filtered cost for the tested payloads as expected, but shows slightly inferior performance to HiLL. With the increase of payload, the performance gap between GMRF and HiLL becomes smaller and GMRF tends to exhibit comparable performance with HiLL for relatively large payload (R$$\geq$$0.3{\rm bpp}). For maxSRMd2, both GMRF and MiPOD show superior performance to HiLL, especially at small payload (R$$\leq$$0.2{\rm bpp}), and the weakness of GMRF compared to MiPOD for relatively large payload (R$$>$$0.3{\rm{bpp}}) is decreasing compared with the results in Table I, which is due to the smoothing operation spreads out the change probabilities.

In short, the smoothing operations indeed boost the security performance of both GMRF and MiPOD against SRM and maxSRMd2, especially for maxSRMd2. On the other hand, the smoothing operation would attenuate the adaptability of GMRF, which is beneficial to the performance against maxSRMd2, and detrimental to the one against SRM. This explains why the performance improvement of GMRF (due to smoothing operation) for SRM is less than the one for maxSRMd2. Finally, although the proposed GMRF (with and without filtering) shows inferior performance to MiPOD for maxSRMd2 at relatively large payload, it exhibits superior security performance consistently to MiPOD for SRM. Considering the fact that, in practice, the precise knowledge of selection channels is generally unavailable to the steganalyzers, so the proposed GMRF is more preferable for practical applications compared with MiPOD.

V Conclusion

At present, the prevailing methodology for adaptive image steganography is based on the framework of minimal distortion embedding, which includes the additive embedding cost for each cover element and the encoding method, typically syndrome-trellis codes (STCs), to minimize the sum of costs. Inspired by the recognition that the security performance of image steganography could be improved by taking advantages of the non-additive model, a Gaussian Markov Random Field (GMRF) with four-element cross neighborhood is proposed to capture the dependences among spatially adjacent pixels, and the problem of secure image steganography is formulated as the minimization of KL-divergence between cover and stego based on sound mathematical principle. The adoption of the proposed GMRF allows to effectively characterize the high-dimensional joint distribution of cover elements and the corresponding KL-divergence in terms of a series of low-dimensional clique structures. With the proposed GMRF, the cover image is tessellated into two disjoint subimages, which are conditionally independent. An alternating iterative optimization scheme is then developed to tackle the issue of efficient embedding while minimizing the total KL-divergence. Finally, the performance of the proposed GMRF is further boosted with smoothing operations on obtained costs. Experiments are carried out to demonstrate the superior performance of the proposed GMRF in terms of secure payload against steganalysis and show that the GMRF outperforms the prior arts, e.g., MiPOD, which is based on the multivariate independent Gaussian model, and has comparable performance with the state-of-the-art HiLL against SRM for tested payloads, where the selection-channel knowledges are unavailable to the steganalyzers and is more preferable for practical applications.

Appendix A KL-divergence and Fisher information matrix for a two-pixel clique

1. KL-divergence

Let ${F_{P}({{\bf X}})}$ and $F_{Q}^{\bm{\beta}}({\bf Y})$ be the p.m.f of two-pixel cover clique ${\bf X}=[X_{1},X_{2}]^{T}$ and stego clique ${\bf Y}=[Y_{1},Y_{2}]^{T}$, respectively. And the associated change probabilities are ${\bm{\beta}}=[{\beta}_{1},{\beta}_{2}]^{T}$, the KL-divergence $D(\bm{\beta})$ between ${\bf X}$ and ${\bf Y}$ is

[TABLE]

By taking the Taylor expansion at ${\bm{\beta}}=[0,0]^{T}$, we have

[TABLE]

Note that $D({\bm{0}})=0$ and $\nabla D({\bm{0}})=0$, for small ${\bm{\beta}}$, $D({\bm{\beta}})$ can then be well approximated with its leading quadratic term

[TABLE]

where ${\nabla^{2}}D({\bm{0}})$ is the second-order partial derivatives defined as

[TABLE]

We then proceed to the computation for each of the components of ${\nabla^{2}}D({\bm{0}})$. Note that ${\left.{F_{Q}^{\bm{\beta}}({\bf{X}}={\bf{x}})}\right|_{{\bf{x}}=[i,j]}}={q_{i,j}}$, according to (10), it is easily verified that

[TABLE]

According to the definition in (A-1), we have

[TABLE]

Note that ${\left.{{F_{P}}({\bf{X}})=F_{Q}^{\bm{\beta}}({\bf{X}})}\right|_{{\bm{\beta}}={\bf{0}}}}$, substituting (A-5) into (A-6) gives

[TABLE]

and similarly we can obtain

[TABLE]

As for the two terms $\displaystyle\frac{{{\partial^{2}}D({\bm{\beta}})}}{{\partial{\beta_{1}}\partial{\beta_{2}}}}$ and $\displaystyle\frac{{{\partial^{2}}D({\bm{\beta}})}}{{\partial{\beta_{2}}\partial{\beta_{1}}}}$, we have

[TABLE]

where

[TABLE]

Similarly, for symmetric embedding, we can verify that (A-10) is equal to [math]. Therefore

[TABLE]

The ${I_{2}}{({\bf{0}})_{1,{\rm{1}}}}$, ${I_{2}}{({\bf{0}})_{1,{\rm{2}}}}$, ${I_{2}}{({\bf{0}})_{2,{\rm{1}}}}$ and ${I_{2}}{({\bf{0}})_{2,{\rm{2}}}}$ in (A-7), (A-8) and (A-11) are the components for $2\times 2$ binary Fisher Information Matrix (FIM) ${{\bf{I}}_{2}}({\bf{0}})$ with respect to $\bf{X}$.

2. Fisher information

Recognize that the KL-divergence for a two-pixel clique $\bf{X}$ can be represented in terms of its FIM, we then further seek to derive ${{\bf{I}}_{2}}({\bf{0}})$ for Gaussian distributed $\bf{X}$. Without loss of generality, we only give the derivation of ${I_{2}}{({\bf{0}})_{1,{\rm{2}}}}$, other components of FIM can be obtained in the same way. By definition,

[TABLE]

By substituting (10) into (A-12) and let $\bm{\beta}\to{\bf{0}}$, we can come to that

[TABLE]

Denote $\Omega_{i,j}^{1}=({p_{i-1,j}}+{p_{i+1,j}}-2{p_{i,j}})$ and ${\Omega_{i,j}^{2}}=({p_{i,j-1}}+{p_{i,j+1}}-2{p_{i,j}})$, we then try to determine (A-13) by taking advantages of the Gaussian distribution $f({\bf X})$ for clique ${\bf X}$. According to (7), $p_{i,j}$ for clique ${\bf{x}}=[x_{1},x_{2}]$ can be formulated as

[TABLE]

where ${f}({x_{1}},{x_{2}})$ is the bivariate Gaussian p.d.f. (see (5)), $i^{\prime}\in(i-0.5,i+0.5)$ and $j^{\prime}\in(j-0.5,j+0.5)$. For fine quantization step $\Delta$, ${f}(i^{\prime}\Delta,j^{\prime}\Delta)\approx{f}(i\Delta,j\Delta)$ , thus we can obtain ${p_{i\pm 1,j}}$ and ${p_{i,j\pm 1}}$ through the Taylor expansion of ${\Delta^{2}}\cdot{f}({x_{1}},{x_{2}})$ at ${x_{1}}=i\Delta$ and ${x_{2}}=j\Delta$ , i.e.,

[TABLE]

[TABLE]

Based on these, $\Omega_{i,j}^{1}$ can then be formulated as

[TABLE]

Therefore, (A-17) can be simplified as

[TABLE]

and similarly, we have

[TABLE]

Substitute (A-13) into (A-12) and note that $q_{i,j}=p_{i,j}$ when ${\bm{\beta}}=\bf{0}$, we can finally obtain

[TABLE]

Similarly, the other two components of ${\bf I}_{2}(\bf 0)$ are

[TABLE]

[TABLE]

Proof for** (A-10):** As an supplement, refer to (A-14), (A-15) and (A-16), we can further obtain ${p_{i\pm 1,j\pm 1}}={F_{\Delta}}((i\pm 1)\Delta,(j\pm 1)\Delta)$ using Taylor expansion of ${\Delta^{2}}\cdot{f}({x_{1}},{x_{2}})$ at ${x_{1}}=i\Delta$ and ${x_{2}}=j\Delta$, simultaneously. For simplicity, we expand it to the second order, and take ${p_{i{\rm{+}}1,j+1}}={F_{\Delta}}((i+1)\Delta,(j+1)\Delta)$ for example, i.e.,

[TABLE]

Similarly, we can get ${p_{i-1,j-1}},{p_{i-1,j{\rm{+}}1}},{p_{i{\rm{+}}1,j-1}}$, then we have

[TABLE]

Next, according to (A-15) and (A-16), we have

[TABLE]

Finally, we substitute (A-24) and (A-25) in (A-10), then we can obtain that the result is [math].

Appendix B KL-divergence for a 4-ary clique tree

For the GMRF model with four-element cross neighborhood as shown in Fig. 1, the four neighboring pixels $\left\{{{X_{1}},{X_{2}},{X_{3}},{X_{4}}}\right\}$ of ${X_{s}}$ are mutually independent when ${X_{s}}$ is given, which constitutes a 4-ary clique tree ${T_{s}}$ centered at ${X_{s}}$. And the joint p.m.f. for ${T_{s}}$ is determined as

[TABLE]

where ${{\cal C}}=\left\{{\{{X_{s}},{X_{t}}\}|t=1,\cdots,4}\right\}$, $F({{\bf{X}}_{c}})$ and $F({X_{s}})$ are the p.m.f. for clique $c=\{{X_{s}},{X_{t}}\}$ and pixel ${X_{s}}$, respectively. For cover clique tree , the embedding modifies it to the stego one with change probabilities ${\bm{\beta}_{{T_{s}}}}=\left\{{{\beta_{X_{s}}},{\beta_{X_{1}}},{\beta_{X_{2}}},{\beta_{X_{3}}},{\beta_{X_{4}}}}\right\}$, by applying (B-1), the KL-divergence between cover and stego clique tree can be written as

[TABLE]

For the first term ${\left\langle 1\right\rangle}$ in (B-2), we have

[TABLE]

For the summation $\sum\limits_{c\in{\cal C}}{\log\displaystyle\frac{{{F_{P}}({{\bf{X}}_{c}})}}{{F_{Q}^{{{\bm{\beta}}_{c}}}{\rm{(}}{{\bf{X}}_{c}}{\rm{)}}}}}$ in (B-3), we take term associated with clique $c=\{{X_{s}},{X_{1}}\}$ for example, other terms can be obtained similarly. We have

[TABLE]

where $\Psi=\left({\sum\limits_{{x_{2}}}{\sum\limits_{{x_{3}}}{\sum\limits_{{x_{4}}}{\displaystyle\frac{{{F_{P}}({X_{s}},{X_{2}}){F_{P}}({X_{s}},{X_{3}}){F_{P}}({X_{s}},{X_{4}})}}{{{{\left({{F_{P}}({X_{s}})}\right)}^{3}}}}}}}}\right)$. It is easily verified that $\sum\limits_{{X_{t}}}{{F_{P}}({X_{s}},{X_{t}})={F_{P}}({X_{s}})}$, then $\Psi=1$, thus (B-4) can be simplified as the KL-divergence of clique $c=\{{X_{s}},{X_{1}}\}$, that is

[TABLE]

On the other hand, the term $\left\langle 2\right\rangle$ in (B-2) can also be simplified as the KL-divergence of pixel ${X_{s}}$ according to [7], i.e.,

[TABLE]

Finally, the KL-divergence for 4-ary clique tree $T_{s}$ can then be written as

[TABLE]

For practical steganography with dynamical clique allocation (see Section III-C), the KL-divergence in (B-7) should be rewritten as

[TABLE]

where $[{\theta_{c}}]=1$ when clique $c$ is included in $T_{s}$ and zero otherwise.

Appendix C Minimization of total KL-divergence between cover and Stego with payload constraint

The minimum distortion embedding with the proposed GMRF model can be formulated as the total KL-divergence minimization of sub-images ${\cal A}$ and ${\cal B}$ with payload constraints. We take the optimization on sub-image ${\cal A}$ for example, by applying the method of Lagrange multipliers, we have

[TABLE]

where $D_{KL}^{{{\cal A}},{T_{s}}}$ is the KL-divergence for $s^{th}$ 4-ary clique tree $T_{s}$ centered at pixel $X_{s}$ in ${\cal A}$, ${\beta_{{X_{s}}}^{{\cal A}}}$ is the change probability for $X_{s}$ which we are going to optimize, and $L/2$ is the payload assigned to ${\cal A}$. Differentiating the objective function with respect to ${\beta_{{X_{s}}}^{{\cal A}}}$ gives

[TABLE]

For term $\left\langle 1\right\rangle$ in (C-2), substituting (B-8) into $\left\langle 1\right\rangle$ gives

[TABLE]

where $\beta_{X_{s\_c}}^{{{\cal B}},{T_{s}}}$ is the change probability for pixel $X_{s\_c}$ corresponding to the clique $c=\{{X_{s}},{X_{s\_c}}\}$ of 4-ary clique tree $T_{s}$, which is located in sub-image ${\cal B}$ and is fixed in the optimization process, ${\bf{I}}_{2,c}^{{T_{s}}}({\bf{0}})$ and $I_{1,X_{s}}^{{T_{s}}}(0)$ are the binary FIM for clique $c$ and FI associated with the $s^{th}$ pixel in ${\cal A}$, respectively. As for the term $\left\langle 2\right\rangle$ in (C-2),we have

[TABLE]

Substituting (C-3) and (C-4) into (C-2), the we can finally obtain

[TABLE]

where ${\Gamma_{{X_{s}}}}=\sum\limits_{c\in{{\cal C}}}{[\theta_{c}^{{T_{s}}}]I_{2,c}^{{T_{s}}}{{({\bf{0}})}_{2,2}}-(}(\sum\limits_{c\in{{\cal C}}}{[\theta_{c}^{{T_{s}}}]})-1)I_{1,X_{s}}^{{T_{s}}}(0)$ and ${\Lambda_{{X_{s}}}}=\sum\limits_{c\in{{\cal C}}}{[\theta_{c}^{{T_{s}}}]I_{2,c}^{{T_{s}}}{{({\bf{0}})}_{1,2}}}\beta_{{X_{s\_c}}}^{{\cal B}}$. And (C-5) could be solved numerically for each $X_{s}$.