Asked by:
For the interested  Fast Poisson Blending for oval regions

Fast Poisson Blending for oval regions
First: If you are looking for an exact solution, or a solution with significant precision at the boundary, dont waste your time with reading this. this is about approximating a "good" result, where "good" means that the visual impression is ok for the viewer of the pictures.
Poisson Blending [https://en.wikipedia.org/wiki/Gradientdomain_image_processing] is one way to seamlessly blend parts of one image into another. Its based on the process of restoring an image from its gradient and some "start" values. If you have an area of a source image plus the gradient of that area, you could replace an equal area in a target image that way, that you take the boundary of the area as fixed start values to restore the colors of the source images gradient in that area of the target image.
[Note: I cannot explain all the math here, its simply too much, so I just briefly show what Methods I use. (What you need to know is some Calculus and Calculus of Variations, some strategies to solve differential equations, some Linear Algebra, a bit of Distributions (Generalized Functions) and Integral Transforms.)]
Mathematically you try to minimize the difference between the gradient of the unknown function (area in the target image) and the gradient of the known function (gradient of the area in the source image), and of course, get the values of the unknown function as a solution. The EulerLagrangeEquation of this Functional is a so called Poisson equation, a partial differential equation.
The whole equation is a Poisson Equation with Dirichlet Boundary conditions (the fixed values of the boundary). Reading [https://www.cs.jhu.edu/~misha/Fall07/Papers/Perez03.pdf] and also this: [http://www.ctralie.com/Teaching/PoissonImageEditing/], explains all the maths and shows a way to retrieve a Matrixvector equation as System to solve the problem. But as you can easily see, the matrix (and vectors) will get huge, since a system for a 200*200 pixels sized area would produce a matrix with 40000 equations and 40000 unknowns.
So its easily imaginable that you can end up with millions of unknowns and equations. The suggested method to solve such a system is an iterative approach like Jacobi or GaussSeidel, and  from my experience  it works quite well up to an area of around 300*300 pixels. But Jacobi and GaussSeidel slow down the more, the larger the system is, and on the other hand it will also take more iterations until a reasonable amount of convergence is reached. So what to do for larger parts?
In terms of image editing we do a deconvolution here. The convolved image is the gradient and we want to restore colors. So why not simply FFT the image and a Laplace Kernel and do the deconvolution in the Fourier domain? There are two reason for not doing so. The first is the Kernel. Transforming a Laplace Kernel (or Gradient seeking Kernels in general) will not produce a "good" shape in the frequency domain, since its truncated very sharply at its edges. So you will get a lot of "ripples" in the fourier domain. (But you can try to use associated functions of fourierpairs directly in the fourier domain!) The second reason is that a simple division  as deconvolution in the fourier domain  will produce very big errors at high frequencies. It can result in an image thats nothing but noise.
But we still can use the FFT here as a Helper for a different type of transform, the discrete sine transform. When inspecting the Matrix of the system, we find, that the Eigenvectors of the Matrix are all sinequations and the corresponding Eigenvalues are of type: 2  (2 * Cos((x+1) * PI / (length+1))) [https://www.youtube.com/watch?v=8XL5ceJZt9Y and https://people.eecs.berkeley.edu/~demmel/cs267/lecture24/lecture24.html]. The Eigenvectors are exactly the imaginary part of the fourier transform if we make our function an odd function, so we can use the FFT for quadratic areas with some pre and post processing to get the correct Sine Transform [https://www5.in.tum.de/lehre/vorlesungen/asc/ss10/vorlesung/dst.pdf]. (If you inspect the Matrix further, you will see that the set of Eigenvectors is complete and all Eigenvalues are real and positive, so the Eigenvectors form an orthonormal Basis of a Vectorspace. Thats the target space for our Transformation). When now dividing by the Eigenvalues and transforming back you will get a result with a lot less high frequency error and its really fast. But, by default, you can use it only for quadratic regions.
So what to do with an oval, or maybe a star shaped region?
The idea is to: a) get the outline of the shape b) project the boundary of the shape to the boundary of the quadratic area and solve by the Fast Sine Transform plus dividing. c) crop the result, so only the initial shape will be left as result.
a) can be done with a floodfill or chaincode implementation.
c) can be don with filling a Graphics Path with a TextureBrush
and b) is just some geometric operations: > get the center point, get the angle of the current (initial boundary) point relative to the center, elongate the line from the center to our current point until it reaches the boundary and write the value of the current coordinate (which is a value from the initial boundary) to the "new" boundary of the quadratic region. Leave all points in between the new boundary and the inner area gradient with a value of zero. Thats the complete setup. The fixed (boundary) conditions are not changed, the inner area is unchanged and the area between the boundary of the quadratic area (that now holds the initial boundary values) and the inner region is zero, so the values wont change when going (solving the system) from the boundary to the inner area. But, when looking left or right from the initial boundary, there might be wrong values, coming from the now rectangular boundary.
The results are quite good, but dont expect too much. this can only be an approximation, since doing so, there will be values in the computation, where there shouldnt. Using the results as input/startVector for some postprocessing Jacobi iterations doesnt work too well, the iterative solver still needs more than half of the amount of iterations, compared to doing it only with a Jacobi solver, when the boundary and the gradient of the upper image are changing much and also when the shape differs a lot from an oval shape. If you use it with eg elliptic shapes and gradients that are relatively flat in the boundary area, you can get very good fft solutions, that will take only a couple of jacobi iterations to a good convergence.
What I'll do next is to *try* to find a way to adjust the projected boundary values by its neighbors values in the initial boundary, so that maybe the results get better. Especially when the boundary itself has some significant/sharp color change in it. The other option would be to modify the gradient values to "work" against producing the errors.
With some post processing (a slight feathering of the edges), the region fits well into the target image. Note that distributing the Boundary values so that there are no zero valued points in the rectanglular boundary can give a slight improvement for images that have an outline of almost the same color, since then it wouldnt add to the high frequencies that much, the most then would go to the low and lowest frequencies. Also, distributing the boundary values can improve results for nonstarshaped regions slightly.
left: pre feathered, FFT, and lower high colors 45, feathered; middle: pre feathered, FFT, feathered; right: pre feathered, FFT, DC at 0.9, feathered
left: pre feathered, Jacobi, 10000 iterations at error 0.001; right: pre feathered, Jacobi, 14000 iterations at error 0.001;
Demo app in VB (PoissonBlendingVB.zip) The App, especially the UI, is at an early stage of development. There might be bugs.
So what's it good for? You'll surely get bored quickly by just blending eyes into hands or sharks into swimmingpools etc. But you could use [an advanced/modified version of] it for blending faint structures into target images. E.G. you could poissonblend a person's head into a target image to get all fine hairstructures and then overlay that with a second image of that person, now well cropped, maybe further processed. Or you could return the blended part as a separate image *for* further processing... Or you could use it to clone colors...
Some more examples:
Edit:
I added a Multigrid method (FMG) to have a relatively fast iterative method for any boundary shapes (not oval, irregular).
Now, since the FFT based method is particularely good for expanding low frequencies for any curved shape and the iterative methods are especially good in expanding the high frequencies, both methods
fit together really good. If you have a shape that doesnt have a rectangular boundary but also doesnt have a convex boundary, you still can use the FFT (Fast Sine Transform) when shifting the boundary to the bounding rectangle, and pass the results as startvectorto the Multigtrid method to get very good results in a reasonable amount of time.
I implemented a method that uses partitioned images (kind of a Finite Tiles Method) for being able to process large images without running out of memory.
> a tiled method with overlapping tiles.
> a tiled method with minimum overlap
continued in next post, for adding more pictures...
Regards,
Thorsten
 Edited by Thorsten Gudera Thursday, February 22, 2018 3:24 PM
General discussion
All replies


... I added the minimum overlap version to the program (standard is 2 pixels, can be changed):
Here's how it works:
As you can see we:
 interpolate on the mesh wire frames
 solve as exact as possible on the Elements (Tiles)
(in contrast to a Finite Element Method)
Regards,
Thorsten
 Edited by Thorsten Gudera Thursday, February 22, 2018 3:22 PM
