Throughout the rest of the paper, we represent images and mathematical models in discrete notations (i.e., vectors and matrices). An image is represented as an M × N matrix, and hence the image domain is denoted by Ω = {1, 2, …, M} × {1, 2, …, N}. We define two inner product spaces: X: = ℝ^M×N and Y: = X × X. Let u ∈ X. For shorthand notation, we define u ≥ 0 if u_i,j ≥ 0 for all (i,j) ∈ Ω. The discrete gradient operator ∇:X→Y is defined by (∇u)_{i, j} = [(_{∇xu)i, j}, (_{∇yu)i, j}], where

and

The space X is equipped with the standard inner product 〈·, ·〉_X, and Euclidean norm ||·||₂. The space Y has the following inner product and norms: for p = (p₁, p₂) ∈ Y and q = (q₁, q₂) ∈ Y,

For brevity, we omit the subscript X or Y in the inner product when its context is clear.

There are two popular discretizations of total variation: the isotropic TV (Rudin et al., 1992) and the anisotropic TV (Choksi et al., 2011), which are defined by

respectively. This work is based on the weighted difference between anisotropic and isotropic TV (AITV) regularization (Lou et al., 2015c), defined by

for a weighting parameter α ∈ [0, 1]. The range of α ensures the non-negativity of the AITV regularization. Note that anisotropic TV is defined as the ℓ₁ norm of the image gradient [(_∇xu)i,j, (_∇yu)i,j] at the pixel location (i, j) ∈ Ω, while isotropic TV is the ℓ₂ norm on the gradient vector. As a result, AITV can be viewed as the ℓ₁ − αℓ₂ regularization on the gradient vector at every pixel, thereby enforcing sparsity individually at each gradient vector.

Poisson noise follows the Poisson distribution with mean and variance η, whose probability mass function is given by

For a clean image g ∈ X, its intensity value at each pixel g_{i, j} serves as the mean and variance for the corresponding noisy observation f ∈ X defined by

To recover the image g from the noisy image f, we find its maximum a posteriori (MAP) estimation u, which maximizes the probability ℙ(u|f). By Bayes' theorem, we have

It further follows from the definition (6) that

Since Poisson noise is i.i.d. pixelwise, we have

The MAP estimate of ℙ(u|f) is equivalent to its negative logarithm, thus leading to the following optimization problem:

The last term −log ℙ(u_i,j) can be regarded as an image prior or a regularization. For example, Le et al. (2007) considered the isotropic total variation as the image prior and proposed a Poisson denoising model

where log is applied pixelwise and 11 is the matrix whose entries are all 1's. The first term in (8) is a concise notation that is commonly used as a fidelity term for Poisson denoising in various imaging applications (Le et al., 2007; Chan et al., 2014; Wen et al., 2016; Chang et al., 2018; Chowdhury et al., 2020a,b).

A Poisson SaT framework (Chan et al., 2014) consists of two steps. Given a noisy grayscale image f ∈ X corrupted by Poisson noise, the first step is the smoothing step that finds a piecewise-smooth solution u^* from the optimization model:

Then in the thresholding step, K − 1 threshold values τ₁ ≤ τ₂ ≤ … ≤ τ_K−1 are appropriately chosen to segment u^* into K regions, where the kth region is given by

with τ0:=infx∈Ωu*(x). The thresholding step is typically performed by k-means clustering.

The Poisson smoothing, lifting, and thresholding (SLaT) framework (Cai et al., 2017) extends the Poisson SaT framework to color images. For a color image f = (f₁, f₂, f₃) ∈ X × X × X, the model (9) is applied to each color channel f_i for i = 1, 2, 3, thus leading to a smoothed color image u*=(u1*,u2*,u3*). An additional lifting step (Luong, 1993) is performed to transform u^* to (u1′,u2′,u3′) in the Lab space (perceived lightness, red-green, and yellow-blue). The channels in Lab space are less correlated than in RGB space, so they may have useful information for segmentation. The RGB image and the Lab image are concatenated to form the multichannel image û: = (u1*,u2*,u3*,u1′,u2′,u3′), followed by the thresholding stage. Generally, k-means clustering yields K centroids c₁, …, c_K as constant vectors, which are used to form the region

for k = 1, …, K such that Ω_k's are disjoint and ⋃k=1KΩk=Ω.

After the thresholding step for both SaT/SLaT, we define a piecewise-constant approximation of the image f by

where c_{k, ℓ} is the ℓth entry of the constant vector c_k and

Recall that d = 1 when f is grayscale, and d = 3 when f is color.

To establish the solution's existence of the proposed model (11), we start with Lemma 1, a discrete version of Poincaré's inequality (Evans, 2010). In addition, we prove Lemma 2 and Proposition 3, leading to the global existence theorem (Theorem 4).

Lemma 1. There exists a constant C > 0 such that

for every u ∈ X and ū:=1MN∑i=1M∑j=1Nui,j.

Proof. We prove it by contradiction. Suppose there exists a sequence {uk}k=1∞ such that

where ūk=1MN∑i=1M∑j=1N(uk)i,j.. For every k, we normalize each element in the sequence by vk=uk-ūk11||uk-ūk11||2. It is straightforward that

By (13), we have

As {vk}k=1∞ is bounded, there exists a convergent subsequence {vkj}j=1∞ such that vkj→v* for v^* ∈ X. It follows from (15) that ||∇v*||2,1=0. Since ker(∇) = {c11:c ∈ ℝ}, then v^* is a constant vector. However, (14) implies that v̄*=0 and ||v*||2=1. This contradiction proves the lemma.               ⃣

Lemma 2. Suppose ||f||_∞ < ∞ and mini,jfi,j>0. There exists a scalar u₀ > 0 such that we have 2(x − f_i,j log x) ≥ x for any x ≥ u₀ and (i, j) ∈ Ω.

Proof. For each (i,j) ∈ Ω, we want to show that there exists u_i,j > 0 such that H(x): = x − 2f_i,j log x ≥ 0 for x ≥ u_i,j. Since H(x) is strictly convex and it attains a global minimum at x = 2f_i,j, it is increasing on the domain x > 2f_i,j. Additionally as x dominates log (x) as x → +∞, there exists u_i,j > 2f_i,j > 0 such that ui,jlogui,j ≥ 2fi,j, which implies that H(u_i,j) = u_i,j − 2f_i,j log u_i,j ≥ 0. As a result, for x ≥ u_i,j > 2f_i,j, we obtain x − 2f_i,j log x = H(x) ≥ H(u_i,j) ≥ 0. Define u0:=maxi,jui,j, and hence we have 2(x − f_i,j log x) ≥ x for x ≥ u₀ ≥ u_i,j, ∀(i, j) ∈ Ω.

Proposition 3. Suppose ker(A) ∩ ker(∇) = {0} and {uk}k=1∞⊂X. If {(Auk,∇uk)}k=1∞ is bounded, then {uk}k=1∞ is bounded.

Proof. Since ker(A) ∩ ker(∇) = {0}, we have A11 ≠ 0. Simple calculations lead to

where the last inequality is due to Lemma 1. The boundedness of {Auk}k=1∞ and {∇uk}k=1∞ implies that {ūk}k=1∞ is also bounded by (16). We apply Lemma 1 to obtain

which thereby proves that {uk}k=1∞ is bounded.

Finally, we adapt the proof in Chan et al. (2014) to establish that F has a global minimizer. Theorem 4. Suppose ||f||_∞ < ∞ and mini,jfi,j>0. If λ > 0, μ ≥ 0, α ∈ [0, 1), and ker(A) ∩ ker(∇) = {0}, then F has a global minimizer.

Proof. It is straightforward that ||∇u||_{2, 1} ≤ ||∇u||₁, thus ||∇u||₁ − α||∇u||_{2, 1} ≥ 0 for α ∈ [0, 1). As a result, we have

Given a scalar f > 0, the function G(x) = x − flog(x) attains its global minimum at x = f. Therefore, we have x − f_i,j log x ≥ f_i,j − f_i,j log f_i,j for all x > 0 and (i,j) ∈ Ω, which leads to a lower bound of F(u), i.e.,

As F(u) is lower bounded by F₀, we can choose a minimizing sequence {uk}k=1∞ and hence F(u_k) has a uniform upper bound, denoted by B₁, i.e., F(u_k) < B₁ for all k ∈ ℕ. It further follows from (17) that

which implies that {|〈Auk-flogAuk,11〉|}k=1∞ is uniformly bounded, i.e., there exists a constant B₂ > 0 such that |〈Au_k − f log Au_k, 11〉| < B₂, ∀k. Using these uniform bounds, we derive that

As α < 1, the sequence {∇uk}k=1∞ is bounded.

To prove the boundedness of {Auk}k=1∞, we introduce the notations of x⁺ = max(x, 0) and x⁻ = −min(x, 0) for any x ∈ ℝ. Then x = x⁺ − x⁻. By Lemma 2, there exists u₀ > 0 such that 2(x − f_i,j log x) ≥ x, ∀x ≥ u₀ and (i,j) ∈ Ω. We observe that

This shows that {Auk}k=1∞ is bounded.

Since both {∇uk}k=1∞ and {Auk}k=1∞ are bounded, then {uk}k=1∞ is bounded due to Proposition 3. Therefore, there exists a subsequence {ukn}n=1∞ that converges to some u^* ∈ X. As F is continuous and thus lower semicontinuous, we have

which means that u^* minimizes F.

To minimize (11), we introduce two auxiliary variables v ∈ X and w = (w_x, w_y) ∈ Y, leading to an equivalent constrained optimization problem:

The corresponding augmented Lagrangian is expressed as

where y ∈ X and z = (z_x, z_y) ∈ Y are Lagrange multipliers and β₁, β₂ are positive parameters. We then apply the alternating direction method of multipliers (ADMM) to minimize (19) that consists of the following steps per iteration k:

where σ > 1.

Remark 1. The scheme presented in (21) slightly differs from the original ADMM (Boyd et al., 2011), the latter of which has σ = 1 in (21f). Having σ > 1 increases the weights of the penalty parameters β_{1, k}, β_{2, k} in each iteration k, thus accelerating the numerical convergence speed of the proposed ADMM algorithm. A similar technique has been used in Cascarano et al. (2021), Gu et al. (2017), Storath and Weinmann (2014), Storath et al. (2014), and You et al. (2019).

All the subproblems (21a)–(21c) have closed-form solutions. In particular, the first-order optimality condition for (21a) is

where Δ = −∇^⊤∇ is the Laplacian operator. If ker(A) ∩ ker(∇) = {0}, then β1,kA⊤A-(μ+β2,k)Δ is positive definite and thereby invertible, which implies that (22) has a unique solution u^k+1. By assuming periodic boundary condition for u, the operators Δ and A^⊤A are block circulant (Wang et al., 2008), and hence (22) can be solved efficiently by the 2D discrete Fourier transform F. Specifically, we have the formula

where F-1 is the inverse discrete Fourier transform, the superscript * denotes complex conjugate, the operation ° is componentwise multiplication, and division is componentwise. By differentiating the objective function of (21b) and setting it to zero, we can get a closed-form solution for v_k+1 given by

where the square root, squaring, and division are performed componentwise. Lastly, the w-subproblem (21c) can be decomposed componentwise as follows:

where the proximal operator for ℓ₁ − αℓ₂ on x ∈ ℝⁿ is given by

The proximal operator for ℓ₁ − αℓ₂ has a closed form solution summarized by Lemma 5.

Lemma 5. (Lou and Yan, 2018) Given x ∈ ℝⁿ, β > 0, and α ∈ [0, 1], the optimal solution to (26) is given by one of the following cases:

where ξ = sign(x)° max (|x| − β, 0).

In summary, we describe the ADMM scheme to solve (11) in Algorithm 2.

We establish the subsequential convergence of ADMM described in Algorithm 2. The global convergence of ADMM (Wang et al., 2019) is inapplicable to our model as the gradient operator ∇ is non-surjective, which will be further investigated in future work. For the sake of brevity, we set β = β₁ = β₂ and denote

In addition, we introduce definitions of subdifferentials (Rockafellar and Wets, 2009), which defines a stationary point of a non-smooth objective function.

Definition 6. For a proper function h:ℝⁿ → ℝ∪{+∞}, define dom(h): = {x ∈ ℝⁿ:h(x) < +∞}.

An important property of the limiting subdifferential is its closedness: for any (x_k, v_k) → (x, v) with v_k ∈ ∂h(x_k), if h(x_k) → h(x), then v ∈ ∂h(x).

Lemma 7. Suppose that ker(A) ∩ ker(∇) = {0} and 0 ≤ α < 1. Let {(uk,vk,wk,yk,zk)}k=1∞ be a sequence generated by Algorithm 2. Then, we have

for some constant ν > 0.

Proof. If ker(A) ∩ ker(∇) = {0}, then β0A⊤A+(β0+μ)∇⊤∇ is positive definite, and hence there exists ν > 0 such that

which implies that Lβk(u,vk,wk,yk,zk) is strongly convex with respect to u with parameter ν. Additionally, Lβk(uk+1,v,wk,yk,zk) is strongly convex with respect to v with parameter β₀ ≤ β_k. It follows from Beck (2017), Theorem 5.25, that we have

As w_k+1 is the optimal solution to (21c), it is straightforward to have

Simple calculations by using (21d)–(21e) lead to

Lastly, we have

Combining (29)–(33) together with the fact that βk=σkβ0 for σ > 1, we obtain