Definition

Traditionally, the image roughness is measured based on the intensity difference between adjacent pixels (40):

where U(µ) represents the image roughness penalty and ψ(x) represents the penalty function. w_j,k represents the weighting factor related to the distance between pixels j and k in neighborhood N_j.

We add this image roughness penalty term to the PWLS-TV method to obtain our proposed PWLS-PR method. The associated mathematical formula for the PWLS-PR method can be expressed as follows:

where α is a regularization parameter that controls the balance between the data fidelity and spatial smoothness. Too much α value reduces noise and makes the reconstructed image smoother but also reduces the resolution, so the choice of α value is very important for the PWLS-PR algorithm.

Inspired by the penalized likelihood method of PET image reconstruction using patch-based regularization (41), we propose using a patch associated with each pixel to calculate the image roughness between neighboring pixels j and k.

In this paper, the patch of a pixel is a square area with the pixel as the center, and the sizes of all the patches in an image are the same. We can understand the meaning of patch as shown in Figure 1. The proposed patch-based roughness function is defined as follows (40):

An example of a sampling window designed for a single pixel (the central red pixel). The white window corresponds to a patch of size 9×9, and the gray window corresponds to a patch of size 3×3.

where f_j(µ) represents the feature vector of the intensity values of all pixels in the patch centered on pixel j. The patch-based distance between pixels j and k is calculated by the following formula:

where j_l represents the lth pixel in the patch associated with pixel j. k_l represents the lth pixel in the patch associated with pixel k. m_l represents the total number of pixels in a patch. r_l represents a positive weighting coefficient equal to the normalized inverse spatial distance between pixel j_l and pixel k_l.

From the existing research, we need to obtain a convex penalty function to guarantee the convergence of the proposed algorithm. Then, the penalty function must meet the following three conditions (50):

The penalty function must be differentiable and symmetric everywhere.

The first-order derivative must be nondecreasing:

When x≥0 and 0 <ω^ψ(0) < +∞, the curve must be nonincreasing:

As discussed in a previous article, examples of possible penalty functions include the following: quadratic functions, the Lange function, the Huber function, and hyperbolic functions (40). When the penalty function is a quadratic function, patch-based regularization is equivalent to pixel-based quadratic regularization. However, the main disadvantage of quadratic regularization is that the discontinuity of the image cannot be taken into account, which may result in excessively smooth edges or fine structures in the reconstructed image.

Based on existing research experience, in this paper, the Lange function is used as the penalty function (40,49):

When |x|<<δ, the function approximates the quadratic function, and when |x|>> δ, the function approaches the absolute function.