maj

@ 0,0 +1,50 @@


# Fast marching inpainting




Following the idea that known image information has to be propagated from the contour of the area to inpaint towards its innermost parts, [Alexander Telea's inpainting algorithm][1] uses Sethian's [fast marching method][2] (FFM) to construct and maintain a list of the pixels forming the contour. The area delimited by this band is progressively shrunk while pixels are processed until none remain to be inpainted.




FFM only helps with the order in which pixels are processed, but does not determine how each pixel is going to be actually inpainted. Telea performs a weighted average of all pixels in the neighborhood of the inpainted pixel. The neighborhood is determined by a radius, which value should be close to the thickness of the area to inpainted. The weight function depends on the following factors:


 the distance between a pixel and it neighbors, ie closers neighbors contribute more;


 the level set distance to the original contour, ie neighbors on the same level set (or iso line) contribute more;


 the collinearity of the vector from a pixel to its neighbors and the FFM direction of propagation. This factor will have the effect of extending isophotes (ie lines) reaching the area to inpaint, by giving more weight to neighbors when they are in the axis going from the inpainting pixel in the direction of propagation of the FFM.




[1]: https://www.rug.nl/research/portal/files/14404904/2004JGraphToolsTelea.pdf


[2]: https://math.berkeley.edu/~sethian/2006/Explanations/fast_marching_explain.html




# Python implementation




Our implementation borrows from several sources, including the [OpenCV C++ implementation][3] and [Telea's implementation][4] itself. As advised in the original paper, we first run a FFM in order to compute distances between pixels outside of the mask and the initial mask contour, before running the main FFM that performs the actual inpainting.




Despite closely following the same algorithm, our Python version is considerably slower than the mentioned implementations. Indeed FFM inpainting is not a vectorized algorithm but rather an iterative one, and therefore doesn't fully benefit from using NumPy. In order to keep the processing time under a reasonable amount, we have chosen to only compute the weighted average previously described, dropping the average gradient that is also mentioned in the article and applied in most implementations. This allows for a x6 speed gain while maintaining "goodenough" results, albeit not as smooth.




[3]: https://github.com/opencv/opencv/blob/master/modules/photo/src/inpaint.cpp


[4]: https://github.com/erich666/jgtcode/tree/master/Volume_09/Number_1/Telea2004/AFMM_Inpainting




# Results




*Click for fullscale image*




 Initial image  Pyheal  OpenCV 


 ::  ::  :: 


 [![][im1_in_thumb]][im1_in]  [![][im1_out_thumb]][im1_out]  [![][im1_cv_thumb]][im1_cv] 




[im1_in]: https://raw.githubusercontent.com/olvb/pyheal/master/samples/tulips_in.png


[im1_in_thumb]: https://raw.githubusercontent.com/olvb/pyheal/master/samples/tulips_in.png


[im1_out]: https://raw.githubusercontent.com/olvb/pyheal/master/samples/tulips_out.png


[im1_out_thumb]: https://raw.githubusercontent.com/olvb/pyheal/master/samples/tulips_out.png


[im1_cv]: https://raw.githubusercontent.com/olvb/pyheal/master/samples/tulips_opencv.png


[im1_cv_thumb]: https://raw.githubusercontent.com/olvb/pyheal/master/samples/tulips_opencv.png




 Initial image  Pyheal  OpenCV 


 ::  ::  :: 


 [![][im2_in_thumb]][im2_in]  [![][im2_out_thumb]][im2_out]  [![][im2_cv_thumb]][im2_cv] 




[im2_in]: https://raw.githubusercontent.com/olvb/pyheal/master/samples/lena_in.png


[im2_in_thumb]: https://raw.githubusercontent.com/olvb/pyheal/master/samples/lena_in.png


[im2_out]: https://raw.githubusercontent.com/olvb/pyheal/master/samples/lena_out.png


[im2_out_thumb]: https://raw.githubusercontent.com/olvb/pyheal/master/samples/lena_out.png


[im2_cv]: https://raw.githubusercontent.com/olvb/pyheal/master/samples/lena_opencv.png


[im2_cv_thumb]: https://raw.githubusercontent.com/olvb/pyheal/master/samples/lena_opencv.png




*Samples images from https://homepages.cae.wisc.edu/~ece533/images/*




The Telea algorithm gives satisfying results for narrow masks. One of its niceties is that it can be directly applied to masks containing noncontiguous regions, without any additional code. When used with larger masks or on textured or patterned images, its halfblurring halfstretching effect will however become apparent.


@ 0,0 +1,22 @@


#!/usr/bin/env python3


import argparse


import pyheal


import imageio




parser = argparse.ArgumentParser()


parser.add_argument('in_path', metavar='input_img', type=str,


help='path to input image')


parser.add_argument('mask_path', metavar='mask_img', type=str,


help='path to mask image')


parser.add_argument('out_path', metavar='ouput_img', type=str,


help='path to output image')


parser.add_argument('r', 'radius', metavar='R', nargs=1, type=int, default=[5],


help='neighborhood radius')




args = parser.parse_args()




img = imageio.imread(args.in_path)


mask_img = imageio.imread(args.mask_path)


mask = mask_img[:, :, 0].astype(bool, copy=False)


pyheal.inpaint(img, mask, args.radius[0])


imageio.imwrite(args.out_path, img)


@ 0,0 +1,27 @@


#!/usr/bin/env python3




import os


import pyheal


import cv2




for file in os.scandir('samples/'):


if not file.name.endswith('_in.png'):


continue




img_name = file.name[:len('_in.png')]




in_path = 'samples/' + img_name + '_in.png'


mask_path = 'samples/' + img_name + '_mask.png'


out_path = 'samples/' + img_name + '_out.png'


cv_path = 'samples/' + img_name + '_opencv.png'




in_img = cv2.imread(in_path)


mask_img = cv2.imread(mask_path)


mask = mask_img[:, :, 0].astype(bool, copy=False)


out_img = in_img.copy()




pyheal.inpaint(out_img, mask, 5)


cv2.imwrite(out_path, out_img)




cv_img = cv2.inpaint(in_img, mask_img[:, :, 0], 5, cv2.INPAINT_TELEA)


cv2.imwrite(cv_path, cv_img)


@ 0,0 +1,280 @@


from math import sqrt as sqrt


import heapq


import numpy as np




# flags


KNOWN = 0


BAND = 1


UNKNOWN = 2


# extremity values


INF = 1e6 # dont use np.inf to avoid inf * 0


EPS = 1e6




# solves a step of the eikonal equation in order to find closest quadrant


def _solve_eikonal(y1, x1, y2, x2, height, width, dists, flags):


# check image frame


if y1 < 0 or y1 >= height or x1 < 0 or x1 >= width:


return INF




if y2 < 0 or y2 >= height or x2 < 0 or x2 >= width:


return INF




flag1 = flags[y1, x1]


flag2 = flags[y2, x2]




# both pixels are known


if flag1 == KNOWN and flag2 == KNOWN:


dist1 = dists[y1, x1]


dist2 = dists[y2, x2]


d = 2.0  (dist1  dist2) ** 2


if d > 0.0:


r = sqrt(d)


s = (dist1 + dist2  r) / 2.0


if s >= dist1 and s >= dist2:


return s


s += r


if s >= dist1 and s >= dist2:


return s


# unsolvable


return INF




# only 1st pixel is known


if flag1 == KNOWN:


dist1 = dists[y1, x1]


return 1.0 + dist1




# only 2d pixel is known


if flag2 == KNOWN:


dist2 = dists[y2, x2]


return 1.0 + dist2




# no pixel is known


return INF




# returns gradient for one pixel, computed on 2 pixel range if possible


def _pixel_gradient(y, x, height, width, vals, flags):


val = vals[y, x]




# compute grad_y


prev_y = y  1


next_y = y + 1


if prev_y < 0 or next_y >= height:


grad_y = INF


else:


flag_prev_y = flags[prev_y, x]


flag_next_y = flags[next_y, x]




if flag_prev_y != UNKNOWN and flag_next_y != UNKNOWN:


grad_y = (vals[next_y, x]  vals[prev_y, x]) / 2.0


elif flag_prev_y != UNKNOWN:


grad_y = val  vals[prev_y, x]


elif flag_next_y != UNKNOWN:


grad_y = vals[next_y, x]  val


else:


grad_y = 0.0




# compute grad_x


prev_x = x  1


next_x = x + 1


if prev_x < 0 or next_x >= width:


grad_x = INF


else:


flag_prev_x = flags[y, prev_x]


flag_next_x = flags[y, next_x]




if flag_prev_x != UNKNOWN and flag_next_x != UNKNOWN:


grad_x = (vals[y, next_x]  vals[y, prev_x]) / 2.0


elif flag_prev_x != UNKNOWN:


grad_x = val  vals[y, prev_x]


elif flag_next_x != UNKNOWN:


grad_x = vals[y, next_x]  val


else:


grad_x = 0.0




return grad_y, grad_x




# compute distances between initial mask contour and pixels outside mask, using FMM (Fast Marching Method)


def _compute_outside_dists(height, width, dists, flags, band, radius):


band = band.copy()


orig_flags = flags


flags = orig_flags.copy()


# swap INSIDE / OUTSIDE


flags[orig_flags == KNOWN] = UNKNOWN


flags[orig_flags == UNKNOWN] = KNOWN




last_dist = 0.0


double_radius = radius * 2


while band:


# reached radius limit, stop FFM


if last_dist >= double_radius:


break




# pop BAND pixel closest to initial mask contour and flag it as KNOWN


_, y, x = heapq.heappop(band)


flags[y, x] = KNOWN




# process immediate neighbors (top/bottom/left/right)


neighbors = [(y  1, x), (y, x  1), (y + 1, x), (y, x + 1)]


for nb_y, nb_x in neighbors:


# skip out of frame


if nb_y < 0 or nb_y >= height or nb_x < 0 or nb_x >= width:


continue




# neighbor already processed, nothing to do


if flags[nb_y, nb_x] != UNKNOWN:


continue




# compute neighbor distance to inital mask contour


last_dist = min([


_solve_eikonal(nb_y  1, nb_x, nb_y, nb_x  1, height, width, dists, flags),


_solve_eikonal(nb_y + 1, nb_x, nb_y, nb_x + 1, height, width, dists, flags),


_solve_eikonal(nb_y  1, nb_x, nb_y, nb_x + 1, height, width, dists, flags),


_solve_eikonal(nb_y + 1, nb_x, nb_y, nb_x  1, height, width, dists, flags)


])


dists[nb_y, nb_x] = last_dist




# add neighbor to narrow band


flags[nb_y, nb_x] = BAND


heapq.heappush(band, (last_dist, nb_y, nb_x))




# distances are opposite to actual FFM propagation direction, fix it


dists *= 1.0




# computes pixels distances to initial mask contour, flags, and narrow band queue


def _init(height, width, mask, radius):


# init all distances to infinity


dists = np.full((height, width), INF, dtype=float)


# status of each pixel, ie KNOWN, BAND or UNKNOWN


flags = mask.astype(int) * UNKNOWN


# narrow band, queue of contour pixels


band = []




mask_y, mask_x = mask.nonzero()


for y, x in zip(mask_y, mask_x):


# look for BAND pixels in neighbors (top/bottom/left/right)


neighbors = [(y  1, x), (y, x  1), (y + 1, x), (y, x + 1)]


for nb_y, nb_x in neighbors:


# neighbor out of frame


if nb_y < 0 or nb_y >= height or nb_x < 0 or nb_x >= width:


continue




# neighbor already flagged as BAND


if flags[nb_y, nb_x] == BAND:


continue




# neighbor out of mask => mask contour


if mask[nb_y, nb_x] == 0:


flags[nb_y, nb_x] = BAND


dists[nb_y, nb_x] = 0.0


heapq.heappush(band, (0.0, nb_y, nb_x))






# compute distance to inital mask contour for KNOWN pixels


# (by inverting mask/flags and running FFM)


_compute_outside_dists(height, width, dists, flags, band, radius)




return dists, flags, band




# returns RGB values for pixel to by inpainted, computed for its neighborhood


def _inpaint_pixel(y, x, img, height, width, dists, flags, radius):


dist = dists[y, x]


# normal to pixel, ie direction of propagation of the FFM


dist_grad_y, dist_grad_x = _pixel_gradient(y, x, height, width, dists, flags)


pixel_sum = np.zeros((3), dtype=float)


weight_sum = 0.0




# iterate on each pixel in neighborhood (nb stands for neighbor)


for nb_y in range(y  radius, y + radius + 1):


# pixel out of frame


if nb_y < 0 or nb_y >= height:


continue




for nb_x in range(x  radius, x + radius + 1):


# pixel out of frame


if nb_x < 0 or nb_x >= width:


continue




# skip unknown pixels (including pixel being inpainted)


if flags[nb_y, nb_x] == UNKNOWN:


continue




# vector from point to neighbor


dir_y = y  nb_y


dir_x = x  nb_x


dir_length_square = dir_y ** 2 + dir_x ** 2


dir_length = sqrt(dir_length_square)


# pixel out of neighborhood


if dir_length > radius:


continue




# compute weight


# neighbor has same direction gradient => contributes more


dir_factor = abs(dir_y * dist_grad_y + dir_x * dist_grad_x)


if dir_factor == 0.0:


dir_factor = EPS




# neighbor has same contour distance => contributes more


nb_dist = dists[nb_y, nb_x]


level_factor = 1.0 / (1.0 + abs(nb_dist  dist))




# neighbor is distant => contributes less


dist_factor = 1.0 / (dir_length * dir_length_square)




weight = abs(dir_factor * dist_factor * level_factor)




pixel_sum[0] += weight * img[nb_y, nb_x, 0]


pixel_sum[1] += weight * img[nb_y, nb_x, 1]


pixel_sum[2] += weight * img[nb_y, nb_x, 2]




weight_sum += weight




return pixel_sum / weight_sum




# main inpainting function


def inpaint(img, mask, radius=5):


if img.shape[0:2] != mask.shape[0:2]:


raise ValueError("Image and mask dimensions do not match")




height, width = img.shape[0:2]


dists, flags, band = _init(height, width, mask, radius)




# find next pixel to inpaint with FFM (Fast Marching Method)


# FFM advances the band of the mask towards its center,


# by sorting the area pixels by their distance to the initial contour


while band:


# pop band pixel closest to initial mask contour


_, y, x = heapq.heappop(band)


# flag it as KNOWN


flags[y, x] = KNOWN




# process his immediate neighbors (top/bottom/left/right)


neighbors = [(y  1, x), (y, x  1), (y + 1, x), (y, x + 1)]


for nb_y, nb_x in neighbors:


# pixel out of frame


if nb_y < 0 or nb_y >= height or nb_x < 0 or nb_x >= width:


continue




# neighbor outside of initial mask or already processed, nothing to do


if flags[nb_y, nb_x] != UNKNOWN:


continue




# compute neighbor distance to inital mask contour


nb_dist = min([


_solve_eikonal(nb_y  1, nb_x, nb_y, nb_x  1, height, width, dists, flags),


_solve_eikonal(nb_y + 1, nb_x, nb_y, nb_x + 1, height, width, dists, flags),


_solve_eikonal(nb_y  1, nb_x, nb_y, nb_x + 1, height, width, dists, flags),


_solve_eikonal(nb_y + 1, nb_x, nb_y, nb_x  1, height, width, dists, flags)


])


dists[nb_y, nb_x] = nb_dist




# inpaint neighbor


pixel_vals = _inpaint_pixel(nb_y, nb_x, img, height, width, dists, flags, radius)




img[nb_y, nb_x, 0] = pixel_vals[0]


img[nb_y, nb_x, 1] = pixel_vals[1]


img[nb_y, nb_x, 2] = pixel_vals[2]




# add neighbor to narrow band


flags[nb_y, nb_x] = BAND


# push neighbor on band


heapq.heappush(band, (nb_dist, nb_y, nb_x))

After Width:  Height:  Size: 544 KiB 
After Width:  Height:  Size: 2.4 KiB 
After Width:  Height:  Size: 508 KiB 
After Width:  Height:  Size: 508 KiB 
After Width:  Height:  Size: 577 KiB 
After Width:  Height:  Size: 2.3 KiB 
After Width:  Height:  Size: 519 KiB 
After Width:  Height:  Size: 519 KiB 
After Width:  Height:  Size: 751 KiB 
After Width:  Height:  Size: 5.4 KiB 
After Width:  Height:  Size: 741 KiB 
After Width:  Height:  Size: 741 KiB 