I implemented Poisson blending in what I think is probably a pretty standard way. First I passed a filter
over the mask to set all pixels without 4 neighbors to zero, so as remove the border of the mask. The remaining
non-zero entries represent the pixels whose values must be calculated, the rest are set equal to the target image.
I created a sparse identity matrix with
a row for each pixel. In the rows that represented unknown values I set the columns that represented the neighbors
of that pixel to -1 and the entry on the diagonal to the number of neighbors. For each color channel, I made
a vector the size of the image and stored the values for all the known pixels as their values in the target image.
The value for an unknown pixel in this vector was set to the difference between that pixel and the average of its
neighbors in the source image. Dividing the sparse matrix by this vector yields the target image with the source
smoothly blended in.
For extra credit, I added a second Poisson function which allows sharp edges from the target to show through
the source. If the difference between the pixel and the average of its neighbors is greater in the target image than
in the source, them the former is used. It didn't seem to make a big difference.