PIL
version of Jack Le's pseudocode for blur with bokeh simulation:
from PIL import ImageFilter
image = Image.open('images/input.png').convert('RGB')
mask = Image.open('images/mask.png').convert('L')
p = 0.4 # power for gamma correction
# r = 3 # in case you want to use Gaussian blur instead of box blur
gammaCorrectedImage = image.point(lambda x: x ** p)
bokeh = gammaCorrectedImage.filter(ImageFilter.BLUR) #ImageFilter.GaussianBlur(radius=r))
bokeh = bokeh.point(lambda x: x ** (1/p)) # To get back the original gamma.
bokeh = Image.composite(image, bokeh, mask)
blurImage = image.filter(ImageFilter.BLUR) #ImageFilter.GaussianBlur(radius=r))
blurImage = Image.composite(image, blurImage, mask)
finalImage = lighter(bokeh, blurImage) #.getdata(), blurImage.getdata())
finalImage.show()
With Box blur kernel and the following input and mask image below
input image

mask image

obtained the following output with the above code (try changing the parameter p
and the blur kernel along with the radius r
to get an output closest to the desired effect):

For example, p=0.5
with Gaussian blur kernel with radius r=7
produces the following output image:
