## 1 Introduction  
Image smoothing is a key technique in image processing and is used in applications such as deblurring [1; 2; 3; 4], detail enhancement [5; 6], tone mapping [7; 8; 9], and so on. Furthermore, it is used in various fields, such as in medical imaging [10; 11; 12], computer graphics [13; 14; 15], remote sensing [16; 17; 18], and character recognition [4; 19; 20]. Filtering-based smoothing methods, e.g., a bilateral filter [6; 21], a guided filter [22; 23; 24; 25], and a nonlocal mean filter [26; 27], have been actively studied for a long time and are often used in practical situations, as they can easily obtain smooth images that roughly maintain the structural gradients of input images with a relatively low computational cost. Optimization-based smoothing methods were recently proposed in [1; 28; 29; 30; 31; 32; 33; 34; 35; 36; 37; 38; 39; 40; 41], and they can flexibly incorporate a priori information into minimization problems. The most standard a priori knowledge of natural images is local smoothness. The popular total variation (TV) [1] is designed as the total magnitude of the vertical and horizontal discrete gradients of an image and promotes this smoothness property on optimization. This method can be employed in more advanced image restoration problems, e.g., reflection removal [42; 43], rain streak removal [44; 45; 46], intrinsic image decomposition [36; 47; 48], image blending [49], and multispectral pansharpening [50; 51; 52; 53], which are complex problems to solve using filtering-based techniques; this is because optimization-based methods can explicitly and quantitatively design observation models by modeling using some norms and specific functions approximately. For example, we can flexibly define not only the commonly used \(\ell_{2}\)-norm but also the \(\ell_{1}\)-norm [54,55] and Huber loss function [49,56,57], which are robust to outliers, and so on, as a data-fidelity term in minimization problems. Among these optimization-based methods, the \(\ell_{0}\) gradient minimization proposed in [31] is known to provide the best approximation of the properties of the local piece-wise flatness of natural images, which is among the most essential prior knowledge in image processing. In this method, the \(\ell_{0}\) gradient, which is the number of nonzero gradients of an image, is minimized together with the \(\ell_{2}\) data fidelity to an input image. Recently, the applications of this method were actively studied, such as algorithm acceleration [58], deblurring [3,4], multi-dimensional data smoothing [58], and layer separation [48,59]. One investigated the \(\ell_{0}\) gradient projection in [37], which allows us to control the degree of smoothing of the output image intuitively. Several extensions of the \(\ell_{0}\) gradient minimization were recently proposed to flexibly separate the structure and texture components in smoothing texture-rich images [39, 40, 41]. However, the \(\ell_{0}\) gradient-based smoothing methods have inherent issues, resulting in undesirable artifacts, such as pseudo-edges and halos in the output images. When an input image contains regions with gradations, the output image often exhibits a piecewise constant appearance, commonly referred to as the staircasing effect. These gradation regions produce strong edges in the output, which are not present in the input image. These issues arise because characteristics related to the intensity and sign of gradients within the observed images are not explicitly considered constraints in optimization problems. Furthermore, these artifacts become more noticeable, particularly in detail enhancement and tone-mapping applications. For example, although these applications are frequently used for diagnosing X-ray images, the staircasing effect introduces a potential risk of misdiagnosis. Similarly, tone mapping is often applied in low-light scenes, such as autonomous driving applications. However, this effect may lead to a reduction in object recognition accuracy.