Modified Poisson compositing technique on skewed grid

1.   Introduction

In image processing, image composing, as an interactive image editing procedure, plays a critical role. Its goal is to create a new image by naturally and seamlessly combining selected image patches from one or more images into one base image. It is commonly used in film making, photo editing, and web design, among other things. Image matting ^[1], image stitching ^[2,3], surface reconstruction ^[4], image completion ^[5], video blending ^[6], and image inpainting ^[7] are examples of other image processing techniques that use a similar approach.

In ^[8], Poisson equation was used to model the image composition process, in which the image gradients were computed with the Laplacian operator. The gradient field inside the cloned region was obtained from the source image, and the Dirichlet boundary conditions were determined by the cloned region where the pixel colour was from the target image. Since then, several works had been conducted in applying image gradients for image processing. Raskar et al. ^[9] used a method that preserves local perceptual qualities while reducing the effects of aliasing, haloing, and ghosting. Sengupta et al. ^[10] suggested a method for solving the Poisson equation in the matte gradient domain to solve the natural image matting problem. Their research showed that the proposed method offered promising results outcomes for a variety of natural images. Wang & Yang ^[2] then developed an image stitching method for combining multiple images with some overlapped images. The similarity between the photos and the visibility of the seam were built in the gradient domain. The authors demonstrated that the proposed strategy was capable of eliminating global inconsistencies and limiting difficulties caused by the illumination effect by working in the gradient domain rather than the pixel intensities of the images. A completion approach based on the gradient domain that used a patch-based technique was presented in ^[5]. The study was divided into two phases: first, they used gradient maps to fill in the missing area using a patch-based filling method; and second, they used the gradient maps to fill in the missing area using a patch-based filling algorithm. By solving the Poisson equation, the image was rebuilt from the gradient maps. A matching criterion based on the gradient and colour was presented to obtain better image completion.

The study in ^[4] for surface reconstruction using an algebraic technique also included image editing based on the gradient domain. They conducted a thorough examination of several integration strategies as well as the difficulty of working with noisy gradient fields. Similar approaches were proposed in ^[11] and ^[12]. The work in ^[13] presented an image composition technique for reducing bleeding artefacts in generated images. In ^[14] the gradient-based approach was applied for video editing. A Fourier implementation to solve the Poisson equation for image processing was reported in ^[15]. The Poisson equation is solved in the Retinex algorithm by Limare et al. ^[16] using Neumann boundary conditions and Fourier implementation. To remove the traces left by camera anonymization used for privacy protection and anti-tracking, Poisson blending was utilized in ^[17]. The suggested strategy efficiently suppresses the visual artefacts produced by camera anonymization while preserving anti-forensic efficiency. Hussain & Kamel ^[18] proposed the image pyramid method for solving the Poisson equation. In their work on image blending, Afifi & Hussain ^[19] suggested the modified Poisson solver. ^[20] and ^[21] proposed several iterative techniques to obtain the solutions of the Poisson equation for image editing purposes.

More recently, Cao et al. ^[22] proposed a self-embedding image watermarking scheme based on reference sharing and the Poisson equation. In their work, they established the relationship between each pixel and its neighborhood in the original image via the Laplacian operator and converted it to compression bits. The approach increased the recovered area from the tampered image by reconstructing the relationship between each compression bit and each reference bit. The work in ^[23] used the Poisson blending technique for cloud removal algorithm. A quad-tree approach was utilized to speed up the proposed gradient-domain compositing method. Experimental photos on a large area of Landsat data with less than 80% cloud coverage yielded encouraging results. Pan et al. ^[24] suggested a method that combines a deep convolutional generative adversarial network and Poisson blending to solve image completion tasks in their work to obtain complete total electron content maps used in global navigation satellite system. In ^[25], a two-stage blending algorithm was presented. In the first stage, a preliminary blended image is synthesized using the Poisson gradient loss, style loss, and content loss. In the second stage, the preliminary blended image is further transformed to have a more similar style to match the target image. It was reported that the approach works well not only for real-world target images but also stylized paintings.

Gradient domain techniques are often required to solve a large sparse linear system, therefore several fast Poisson solvers ^{[21,26,27,28]} were proposed in the past. Despite the advantages of the previous Poisson solvers, the computation becomes time-consuming when applied to large images, where gradient field was obtained on a fine grid. In ^[29], a reduced-order finite difference method based on proper orthogonal decomposition technique was applied.

In this paper, we propose a novel fast Poisson solver that utilized red-black ordering strategy with two different relaxation factors and two accelerated parameters implemented on skewed grid. With red-black ordering, each pixel in the computed region is labeled as red or black color node. Since red and black nodes are computed independently, different relaxation factors can be applied to the two groups of nodes, thus reduces the number of iterations required during the compositing process. The additional relaxation factors and accelerated parameters provide wider choice of values to further speed up the execution time. By implementing the compositing process on a skewed grid, the total computations can be drastically reduced since only half of the pixels inside the compositing region are involved during the iteration procedure. In the previous studies, the iterative algorithms implemented on skewed grid were successfully applied to solve various problems ^[30,31,32].

The main contributions of this article are as follows:

● We propose a new image compositing method implemented on a skewed grid, in which we develop a new algorithm using red-black ordering strategy.

● We execute image compositing process from gradient field obtained by solving the Poisson equation using iterative solvers that employ additional relaxation factor and accelerated parameters for wider parameter tuning options to improve the overall computational time.

● We conduct a series of experiments, and the results demonstrate the superiority of our image compositing method.

The remainder of this paper is organized as follows. Section 2 establishes a Poisson image compositing approach that employ red-black ordering strategy on regular and skewed grids. Experimental evaluations for comparing the proposed solution on skewed grid against solutions on regular grid are presented in Section 3. Discussion on the findings of the suggested algorithms and their limitations are provided in Section 4. Finally, Section 5 presents the conclusions and an outlook on possible future work.

2.   Materials and methods

2.1.   Poisson image compositing

Directly copying a foreground object from source image and pasting it onto a target image can produce big intensity changes at the boundary, which creates obvious artifacts to human eyes. Therefore, the motivation of Poisson image compositing is to smooth the abrupt intensity change in the compositing boundary in order to reduce artifacts. In Figure 1, the left image shows the abrupt intensity change in the composite boundary in the copy-and-paste image and a smooth boundary transition using Poisson compositing ^[8].

Based on the original work in ^[8], Poisson image compositing is formulated as an image interpolation problem using a guidance vector field.

where $\nabla = \left[ \frac{\partial.}{\partial x}, \frac{\partial.}{\partial y} \right]$ is the gradient operator, $f$ is the function of the compositing image, $f^{\ast}$ is the function of the target image, v is the vector field, $\Omega$ is the compositing region and $\partial\Omega$ is the boundary of the compositing region. The solution is the unique solution of the following Poisson equation with Dirichlet boundary conditions:

where $\text{div}{\bf v} = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y}$ is the divergence of ${{\bf v}} = (u, v)$. This is the fundamental machinery of Poisson editing of color images: three Poisson equations of the form (2.1) are solved independently in the three color channels of the chosen color space. All the results reported in this paper were obtained in the RGB color space.

2.1.1.   Iterative solver

This section describes the iterative solver to solve the Poisson equation (2.2). For discrete images, the problem can be discretized naturally using the underlying discrete pixel grid. The finite difference discretization of Eq (2.2) using the five-point Laplacian operator will result in the following approximation equation:

The iterative scheme of Eq (2.3) is given as follows ^[33]:

From Eq (2.4), the Successive Overrelaxation (SOR) ^[33] scheme that utilizes a relaxation factor is obtained and written as follows:

In ^[34], the Accelerated Overrelaxation (AOR) scheme was developed and described as follows:

where an accelerated parameter $\rho$ is added. An extension to the SOR and AOR is the Two-Parameter Overrelaxation (TOR) iterative method ^[35] which employs two accelerated parameters $\rho$ and $\sigma$, to provide more tuning choices during the iteration process. The finite-difference approximation equation for the TOR method is given as follows:

The optimal values for the three parameters ($\omega$, $\rho$ and $\sigma$) are in the range between 1 and 2.

2.2.   Modified iterative algorithm using red-black ordering

As shown in Figure 2, in contrast to the standard ordering in which all nodes are in one color, the modified iterative algorithm employs red-black ordering approach that sets the nodes in red and black colors ^[36]. By applying the red-black ordering, the computations of red and black nodes involve two phases, in which the red and black nodes are computed alternately. The procedure begins by computing red nodes where ($i+j$) is even. In the second phase, black nodes, where ($i+j$) is odd, are computed. As shown in Figure 3, the computations of red and black nodes $u^{(k)}$ are based on the updated values $u^{(k+1)}$ of their four neighbouring opposite color nodes. The iteration procedures of these two phases continue until the stopping criterion is satisfied.

Based on the iteration scheme given in Eq (2.5), the Modified Successive Overrelaxation (MSOR) algorithm ^[37] can be written as

and

where $\alpha$ and $\beta$ are the relaxation factors for red and black nodes respectively, and both values are chosen such that $0 < \alpha, \beta < 2$, where the optimal values are in the range between $1$ and $2$.

The Modified Accelerated Overrelaxation (MAOR) algorithm is based on Eq (2.6), where Eq (2.8) is used for the iteration of its red nodes and the iterative scheme for black nodes is written as

where $0 < \alpha, \beta, \rho < 2$.

Based on Eq (2.7), the Modified Two-Parameter Overrelaxation (MTOR) algorithm employs two relaxation factors and two accelerated parameters. Similar to the MSOR and MAOR, the iteration for its red nodes uses Eq (2.8). Meanwhile, the iterative scheme of black nodes for the MTOR is written as

where all values of relaxation factors ($\alpha$ and $\beta$) and accelerated parameters ($\rho$ and $\sigma$) are between $0$ and $2$.

2.3.   Iterative algorithm on skewed grid

Apart from standard iterative algorithm on regular grid, approximation equation for problem (2.2) can be implemented on skewed grid ^[38] (see Figure 4), in which the approximation equation can be written as

and its iterative scheme is given as

Figure 5 shows the computational molecules of red and black nodes on a skewed grid.

2.3.1.   Skewed modified algorithms

The SOR ^[38], AOR, and TOR can be implemented on skewed grid using difference equation (2.12). Therefore, their corresponding iterative algorithms on skewed grid can be obtained as

and

respectively. By utilizing red-black ordering as illustrated in Figure 6, the skewed modified variants namely Skewed MSOR (SkMSOR), Skewed MAOR (SkMAOR), and Skewed MTOR (SkMTOR) can be obtained. The iterative scheme for these three algorithms on red nodes is written as

Accordingly, the iterative schemes of black nodes for the respective SkMSOR, SkMAOR and SkMTOR are given as

and

2.4.   Image compositing algorithm

The algorithm is only applied to the pixels involved in the compositing process. All other pixels outside the compositing region are ignored during the computation process, thus speed up the overall execution time. By employing the red-black ordering approach, the two different relaxation factors $\alpha$ and $\beta$ can be tuned and thus providing additional options to improve the execution time. Note that if $\alpha = \beta$, the original non-modified variants are obtained. With TOR variants, two accelerated parameters $\rho$ and $\sigma$ are available for tuning, thus provide more choice during the composition process. If $\rho = \sigma$, the AOR variants are obtained, whereas if $\rho = \sigma = 1$, we shall obtain the standard SOR variants. Algorithm 1 describes the iterative procedure for red-black ordering approach, in which the computations for red and black nodes are executed alternately until the maximum tolerance error is obtained. The $\mathtt{ComputePixels}$ function utilizes the appropriate iterative algorithm described in the previous section. For computation on skewed grid, the remaining white nodes are computed after the convergence is achieved.

2.5.   Computational complexity analysis

This section presents the analysis on the total computing costs for the considered algorithms. All algorithms were tested using the same machine running Xubuntu 18.04 on Intel i5 3470K @ 2.6GHz with 16GB memory. By assuming that the number of pixels involve in the calculation is $N$, the total arithmetic operations per iteration can be estimated for each algorithm. On skewed grid, only half of the total pixels, $M = \frac{N}{2}$, are involved during the computations. Meanwhile, for the modified variants that employ red-black ordering approach, the total pixels involved during the iterations of red and black nodes are $M$ on regular grid and $m = \frac{N}{4}$ on skewed grid. We assume that the execution time for the arithmetic addition and multiplication operations require the same computing clock cycle.

For computations on a regular grid, the arithmetic operations required per iteration, for the standard algorithms (SOR, AOR, and TOR) are $7N$, $12N$, and $13N$, respectively. Their corresponding modified algorithms (MSOR, MAOR, and MTOR) require $7M$ for red nodes and $7M$, $16M$, $17M$ for black nodes, respectively. On skewed grid, the respective arithmetic operations per iteration are $7M$, $12M$, and $13M$. The skewed modified algorithms (SkMSOR, SkMAOR, and SkMTOR) require $7m$ arithmetic operations for red nodes. Their respective arithmetic operations for black nodes are $7m$, $16m$, and $17m$. The total computing cost for the considered iterative algorithms are summarized in Table 1.

2.6.   Image quality measurements

The image quality measurement process is used to compare the similarity between the final output images. Based on the statistical method Analysis of Variance (ANOVA) ^[35], seven metrics are used namely Mean Square Error (MSE), Structural Similarity Index (SSIM) ^[39], Structural Content (SC), Normalized Cross-Correlation (NCC), Normalized Absolute Error (NAE), Average Difference (AD) and Maximum Difference (MD) ^[40].

The simple MSE can be written as

where $A$ and $B$ represent the pixel values of reference and target images, and ($m \times n$) is the size of the image. MSE is used to measures the difference between pixel values in $A$ and $B$, in which a smaller value means higher similarity. The ideal MSE value of $0$ is obtained when the two images $A$ and $B$ are identical. Another similarity test that can be used to compare the two images is SSIM ^[41]. If the obtained SSIM is equal to $1$, it means the two images identical.

Similarly, SC and NCC values that are equal to $1$ can be used to indicate that the two images are identical. SC value can be obtained using the following formula:

NCC measurement compares the two images and is given as follows:

In contrast, NAE and AD values that are close to $0$ means the two images are identical. Both values can be obtained as follows:

and

where $A$ and $B$ are the reference and target images, respectively.

Lastly, the MD value that is used to obtain the maximum difference of pixel value between the two images $A$ and $B$. MD is given the following formula:

If the obtained MD value is very small, it means the similarity between the two images is higher. All these seven metrics are used to compare the similarity of the output images generated by the six tested methods.

3.   Results

In this work, three sets of images are used to measure the performance of the considered methods. The six tested methods are measured in terms of the number of iterations, execution time, and quality of images produced. All variants require optimal relaxation factor and accelerated parameter values. By conducting a trial and error procedure, initial results show that all parameter values should be in the range between 1.5 and 1.9. Based on these results, values between 1.60 and 1.74 are chosen, as shown in Table 2.

The quality of the images is compared using the similarity index as described in the previous section. Table 3 shows the number of pixels inside the image masks that were applied for the three image sets. Figure 7 shows all the image sets that comprise source, target, mask, and initial images. In Figure 8, the gradual improvement of the compositing process is illustrated. For sky and ballon images, all regular variants (MSOR, MAOR, and MTOR) take at least 600 iterations to converge. The skewed variants (SkMSOR, SkMAOR, and SkMTOR) reduce the required iterations by approximately 50%.

Table 2 shows the performance of the six tested methods. The proposed SkMTOR method gives the lowest number of iterations, thus consequently provides the best execution time. All modified skewed variants, namely SkMSOR, SkMAOR, and SkMTOR, perform faster than their corresponding regular variants of MSOR, MAOR, and MTOR. It is clearly shown from the results that the SkMTOR outperforms the other methods. Also, since all methods applied red-black ordering, the proposed methods can be implemented for parallel execution.

From Table 4, the MAOR and SkMAOR are better than the respective MSOR and SkMSOR by decreasing the iterations by 8%. The MTOR and SkMTOR slightly decrease the number of iterations by approximately 5% compared to both MAOR and SkMAOR, respectively. Both MTOR and SkMTOR are also clearly superior to the the corresponding MSOR and SkMSOR by approximately 12%. All skewed variants (SkMSOR, SkMAOR, and SkMTOR) outperforms their corresponding regular variants (MSOR, MAOR, and MTOR) by reducing the iterations by approximately 46%. In terms of CPU time, both MAOR and SkMAOR reduce the computational execution by 7% compared to the MSOR and SkMSOR, respectively. The MTOR improve the execution time by 5% compared to the MAOR, whereas the SkMTOR is slightly faster than the SkMAOR by approximately 4%. The MTOR is clearly faster than the MSOR by approximately 12%, while SkMTOR speed up the execution time by 10% compared to the SkMSOR. All skewed variants are clearly faster than their corresponding regular variants by drastically decreasing the execution time by half at approximately 51 to 52%.

As shown in Figure 9, the generated images of grass and bird (top), and lake and crocodile (bottom) obtained from compositing process are all visually identical. No noticeable differences can be seen from all images generated by the six methods. The results of the seven metrics measurement are given in Table 5.

4.   Discussion

The findings of the algorithms suggest that the red-black technique is a viable alternative to existing approaches to tackle image composition problems. The red-black technique is well suited to parallel implementation since the computation for red and black nodes can be done independently, like on a chessboard.

Because only half of the total nodes are engaged in the iteration process, the implementation on a skewed grid minimises computational complexity. As a result, this method dramatically reduces the number of iterations and speeds up the convergence rate.

The MSE, NAE, and AD values are all close to 0 based on the obtained similarity measurement. Meanwhile, the SSIM and NCC are both close to 1, implying that the images generated are nearly identical. The highest difference (MD) between benchmark and generated image pixel intensity values is very low (less than 10). All of these similarity measurement results show that all of the algorithms examined yield identical images.

The proposed composing technique, however, requires manual fine tuning, in which optimal parameter values are obtained through a trial and error procedure. Moreover, the color of the object from the source image becomes darker or brighter, depending upon the background color of the target image.

5.   Conclusions

For the six iterative approaches considered in this study, the red-black ordering implementation is employed to solve the image compositing problem. According to the findings of the experiments, the SkMTOR approach gives the fastest rate of convergence. All skewed variants outperformed their regular versions. The red-black technique, which facilitates parallelism, was not used in most prior studies. The generated images were compared using seven image similarity indexes. The image similarity tests revealed that the suggested methods produced almost identical images with no noticeable differences. More powerful point iterative approaches, such as quarter-sweep iteration, can be used to enhance the work in the future. Block-based iterative approaches would also be studied.

Acknowledgments

The authors would like to thank the Universiti Malaysia Sabah for financially supporting this research under grant SGA0117-2019.

Conflict of interest

The authors declare no conflict of interest.