diff --git a/Watermarking .ipynb b/Watermarking .ipynb index f38a2bd..2f84991 100644 --- a/Watermarking .ipynb +++ b/Watermarking .ipynb @@ -373,21 +373,21 @@ ], "metadata": { "kernelspec": { - "display_name": "Python 2", + "display_name": "Python 3", "language": "python", - "name": "python2" + "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", - "version": 2 + "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", - "pygments_lexer": "ipython2", - "version": "2.7.12" + "pygments_lexer": "ipython3", + "version": "3.8.2" } }, "nbformat": 4, diff --git a/final/fotolia_137840668.jpg b/final/fotolia_137840668.jpg deleted file mode 100644 index 6725ed8..0000000 Binary files a/final/fotolia_137840668.jpg and /dev/null differ diff --git a/final/fotolia_168667147.jpg b/final/fotolia_168667147.jpg deleted file mode 100644 index ef5481d..0000000 Binary files a/final/fotolia_168667147.jpg and /dev/null differ diff --git a/final/fotolia_168667186.jpg b/final/fotolia_168667186.jpg deleted file mode 100644 index b501fe9..0000000 Binary files a/final/fotolia_168667186.jpg and /dev/null differ diff --git a/final/fotolia_168667261.jpg b/final/fotolia_168667261.jpg deleted file mode 100644 index 9afa051..0000000 Binary files a/final/fotolia_168667261.jpg and /dev/null differ diff --git a/final/fotolia_168667468.jpg b/final/fotolia_168667468.jpg deleted file mode 100644 index 036b7a2..0000000 Binary files a/final/fotolia_168667468.jpg and /dev/null differ diff --git a/final/fotolia_168667490.jpg b/final/fotolia_168667490.jpg deleted file mode 100644 index 3d2b4ee..0000000 Binary files a/final/fotolia_168667490.jpg and /dev/null differ diff --git a/final/fotolia_168668046.jpg b/final/fotolia_168668046.jpg deleted file mode 100644 index 90d6eda..0000000 Binary files a/final/fotolia_168668046.jpg and /dev/null differ diff --git a/final/fotolia_168668148.jpg b/final/fotolia_168668148.jpg deleted file mode 100644 index 1721a72..0000000 Binary files a/final/fotolia_168668148.jpg and /dev/null differ diff --git a/final/fotolia_168668150.jpg b/final/fotolia_168668150.jpg deleted file mode 100644 index 4492575..0000000 Binary files a/final/fotolia_168668150.jpg and /dev/null differ diff --git a/final/fotolia_168668190.jpg b/final/fotolia_168668190.jpg deleted file mode 100644 index e384361..0000000 Binary files a/final/fotolia_168668190.jpg and /dev/null differ diff --git a/final/fotolia_75353029.jpg b/final/fotolia_75353029.jpg deleted file mode 100644 index 6524531..0000000 Binary files a/final/fotolia_75353029.jpg and /dev/null differ diff --git a/main.py b/main.py index f678104..7d47b98 100644 --- a/main.py +++ b/main.py @@ -4,44 +4,50 @@ # est = poisson_reconstruct(gx, gy, np.zeros(gx.shape)[:,:,0]) cropped_gx, cropped_gy = crop_watermark(gx, gy) -W_m = poisson_reconstruct(cropped_gx, cropped_gy) # random photo -img = cv2.imread('images/fotolia_processed/fotolia_137840645.jpg') +img = cv2.imread('images/fotolia_processed/fotolia_137840668.jpg') im, start, end = watermark_detector(img, cropped_gx, cropped_gy) +# # Save result of watermark detection +# plt.subplot(1, 2, 1) +# plt.imshow(img) +# plt.subplot(1, 2, 2) # plt.imshow(im) -# plt.show() +# plt.savefig('images/results/watermark_detect_result.png') # We are done with watermark estimation + # W_m is the cropped watermark +W_m = poisson_reconstruct(cropped_gx, cropped_gy) + +# Get the number of images in the folder num_images = len(gxlist) -J, img_paths = get_cropped_images('images/fotolia_processed', num_images, start, end, cropped_gx.shape) +J, img_paths = get_cropped_images( + 'images/fotolia_processed', num_images, start, end, cropped_gx.shape) + # get a random subset of J -idx = [389, 144, 147, 468, 423, 92, 3, 354, 196, 53, 470, 445, 314, 349, 105, 366, 56, 168, 351, 15, 465, 368, 90, 96, 202, 54, 295, 137, 17, 79, 214, 413, 454, 305, 187, 4, 458, 330, 290, 73, 220, 118, 125, 180, 247, 243, 257, 194, 117, 320, 104, 252, 87, 95, 228, 324, 271, 398, 334, 148, 425, 190, 78, 151, 34, 310, 122, 376, 102, 260] -idx = idx[:25] -# Wm = (255*PlotImage(W_m)) -Wm = W_m - W_m.min() # get threshold of W_m for alpha matte estimate -alph_est = estimate_normalized_alpha(J, Wm) +Wm = W_m - W_m.min() + +# estimate alpha matte +alph_est = estimate_normalized_alpha(J, Wm, num_images=len(J)) alph = np.stack([alph_est, alph_est, alph_est], axis=2) C, est_Ik = estimate_blend_factor(J, Wm, alph) alpha = alph.copy() -for i in xrange(3): - alpha[:,:,i] = C[i]*alpha[:,:,i] +for i in range(3): + alpha[:, :, i] = C[i]*alpha[:, :, i] Wm = Wm + alpha*est_Ik W = Wm.copy() -for i in xrange(3): - W[:,:,i]/=C[i] +for i in range(3): + W[:, :, i] /= C[i] Jt = J[:25] + # now we have the values of alpha, Wm, J # Solve for all images Wk, Ik, W, alpha1 = solve_images(Jt, W_m, alpha, W) -# W_m_threshold = (255*PlotImage(np.average(W_m, axis=2))).astype(np.uint8) -# ret, thr = cv2.threshold(W_m_threshold, 127, 255, cv2.THRESH_BINARY) - diff --git a/src/__init__.py b/src/__init__.py index b047e59..e5ce21b 100644 --- a/src/__init__.py +++ b/src/__init__.py @@ -1,4 +1,4 @@ -from estimate_watermark import * -from preprocess import * -from image_crawler import * -from watermark_reconstruct import * +from .estimate_watermark import * +from .preprocess import * +from .image_crawler import * +from .watermark_reconstruct import * diff --git a/src/estimate_watermark.py b/src/estimate_watermark.py index 3ee5425..fa3557e 100644 --- a/src/estimate_watermark.py +++ b/src/estimate_watermark.py @@ -1,182 +1,196 @@ -import sys, os +import sys +import os import cv2 import numpy as np import warnings from matplotlib import pyplot as plt import math import numpy -import scipy, scipy.fftpack +import scipy +import scipy.fftpack # Variables KERNEL_SIZE = 3 + def estimate_watermark(foldername): - """ - Given a folder, estimate the watermark (grad(W) = median(grad(J))) - Also, give the list of gradients, so that further processing can be done on it - """ - if not os.path.exists(foldername): - warnings.warn("Folder does not exist.", UserWarning) - return None - - images = [] - for r, dirs, files in os.walk(foldername): - # Get all the images - for file in files: - img = cv2.imread(os.sep.join([r, file])) - if img is not None: - images.append(img) - else: - print("%s not found."%(file)) - - # Compute gradients - print("Computing gradients.") - gradx = map(lambda x: cv2.Sobel(x, cv2.CV_64F, 1, 0, ksize=KERNEL_SIZE), images) - grady = map(lambda x: cv2.Sobel(x, cv2.CV_64F, 0, 1, ksize=KERNEL_SIZE), images) - - # Compute median of grads - print("Computing median gradients.") - Wm_x = np.median(np.array(gradx), axis=0) - Wm_y = np.median(np.array(grady), axis=0) - - return (Wm_x, Wm_y, gradx, grady) + """ + Given a folder, estimate the watermark (grad(W) = median(grad(J))) follow Eq.4 + Also, give the list of gradients, so that further processing can be done on it + """ + if not os.path.exists(foldername): + warnings.warn("Folder does not exist.", UserWarning) + return None + + images = [] + for r, dirs, files in os.walk(foldername): + # Get all the images + for file in files: + img = cv2.imread(os.sep.join([r, file])) + if img is not None: + images.append(img) + else: + print("%s not found." % (file)) + + # Compute gradients + print("Computing gradients.") + gradx = list(map(lambda x: cv2.Sobel( + x, cv2.CV_64F, 1, 0, ksize=KERNEL_SIZE), images)) + grady = list(map(lambda x: cv2.Sobel( + x, cv2.CV_64F, 0, 1, ksize=KERNEL_SIZE), images)) + + # Compute median of grads + print("Computing median gradients.") + Wm_x = np.median(np.array(gradx), axis=0) + Wm_y = np.median(np.array(grady), axis=0) + + return (Wm_x, Wm_y, gradx, grady) def PlotImage(image): - """ - PlotImage: Give a normalized image matrix which can be used with implot, etc. - Maps to [0, 1] - """ - im = image.astype(float) - return (im - np.min(im))/(np.max(im) - np.min(im)) + """ + PlotImage: Give a normalized image matrix which can be used with implot, etc. + Maps to [0, 1] + """ + im = image.astype(float) + return (im - np.min(im))/(np.max(im) - np.min(im)) def poisson_reconstruct2(gradx, grady, boundarysrc): - # Thanks to Dr. Ramesh Raskar for providing the original matlab code from which this is derived - # Dr. Raskar's version is available here: http://web.media.mit.edu/~raskar/photo/code.pdf - - # Laplacian - gyy = grady[1:,:-1] - grady[:-1,:-1] - gxx = gradx[:-1,1:] - gradx[:-1,:-1] - f = numpy.zeros(boundarysrc.shape) - f[:-1,1:] += gxx - f[1:,:-1] += gyy - - # Boundary image - boundary = boundarysrc.copy() - boundary[1:-1,1:-1] = 0; - - # Subtract boundary contribution - f_bp = -4*boundary[1:-1,1:-1] + boundary[1:-1,2:] + boundary[1:-1,0:-2] + boundary[2:,1:-1] + boundary[0:-2,1:-1] - f = f[1:-1,1:-1] - f_bp - - # Discrete Sine Transform - tt = scipy.fftpack.dst(f, norm='ortho') - fsin = scipy.fftpack.dst(tt.T, norm='ortho').T - - # Eigenvalues - (x,y) = numpy.meshgrid(range(1,f.shape[1]+1), range(1,f.shape[0]+1), copy=True) - denom = (2*numpy.cos(math.pi*x/(f.shape[1]+2))-2) + (2*numpy.cos(math.pi*y/(f.shape[0]+2)) - 2) - - f = fsin/denom - - # Inverse Discrete Sine Transform - tt = scipy.fftpack.idst(f, norm='ortho') - img_tt = scipy.fftpack.idst(tt.T, norm='ortho').T - - # New center + old boundary - result = boundary - result[1:-1,1:-1] = img_tt - - return result - - -def poisson_reconstruct(gradx, grady, kernel_size=KERNEL_SIZE, num_iters=100, h=0.1, - boundary_image=None, boundary_zero=True): - """ - Iterative algorithm for Poisson reconstruction. - Given the gradx and grady values, find laplacian, and solve for image - Also return the squared difference of every step. - h = convergence rate - """ - fxx = cv2.Sobel(gradx, cv2.CV_64F, 1, 0, ksize=kernel_size) - fyy = cv2.Sobel(grady, cv2.CV_64F, 0, 1, ksize=kernel_size) - laplacian = fxx + fyy - m,n,p = laplacian.shape - - if boundary_zero == True: - est = np.zeros(laplacian.shape) - else: - assert(boundary_image is not None) - assert(boundary_image.shape == laplacian.shape) - est = boundary_image.copy() - - est[1:-1, 1:-1, :] = np.random.random((m-2, n-2, p)) - loss = [] - - for i in xrange(num_iters): - old_est = est.copy() - est[1:-1, 1:-1, :] = 0.25*(est[0:-2, 1:-1, :] + est[1:-1, 0:-2, :] + est[2:, 1:-1, :] + est[1:-1, 2:, :] - h*h*laplacian[1:-1, 1:-1, :]) - error = np.sum(np.square(est-old_est)) - loss.append(error) - - return (est) + # Thanks to Dr. Ramesh Raskar for providing the original matlab code from which this is derived + # Dr. Raskar's version is available here: http://web.media.mit.edu/~raskar/photo/code.pdf + + # Laplacian + gyy = grady[1:, :-1] - grady[:-1, :-1] + gxx = gradx[:-1, 1:] - gradx[:-1, :-1] + f = numpy.zeros(boundarysrc.shape) + f[:-1, 1:] += gxx + f[1:, :-1] += gyy + + # Boundary image + boundary = boundarysrc.copy() + boundary[1:-1, 1:-1] = 0 + + # Subtract boundary contribution + f_bp = -4*boundary[1:-1, 1:-1] + boundary[1:-1, 2:] + \ + boundary[1:-1, 0:-2] + boundary[2:, 1:-1] + boundary[0:-2, 1:-1] + f = f[1:-1, 1:-1] - f_bp + + # Discrete Sine Transform + tt = scipy.fftpack.dst(f, norm='ortho') + fsin = scipy.fftpack.dst(tt.T, norm='ortho').T + + # Eigenvalues + (x, y) = numpy.meshgrid( + range(1, f.shape[1]+1), range(1, f.shape[0]+1), copy=True) + denom = (2*numpy.cos(math.pi*x/(f.shape[1]+2))-2) + \ + (2*numpy.cos(math.pi*y/(f.shape[0]+2)) - 2) + + f = fsin/denom + + # Inverse Discrete Sine Transform + tt = scipy.fftpack.idst(f, norm='ortho') + img_tt = scipy.fftpack.idst(tt.T, norm='ortho').T + + # New center + old boundary + result = boundary + result[1:-1, 1:-1] = img_tt + + return result + + +def poisson_reconstruct(gradx, grady, kernel_size=KERNEL_SIZE, num_iters=100, h=0.1, + boundary_image=None, boundary_zero=True): + """ + Iterative algorithm for Poisson reconstruction. + Given the gradx and grady values, find laplacian, and solve for image + Also return the squared difference of every step. + h = convergence rate + """ + fxx = cv2.Sobel(gradx, cv2.CV_64F, 1, 0, ksize=kernel_size) + fyy = cv2.Sobel(grady, cv2.CV_64F, 0, 1, ksize=kernel_size) + laplacian = fxx + fyy + m, n, p = laplacian.shape + + if boundary_zero == True: + est = np.zeros(laplacian.shape) + else: + assert(boundary_image is not None) + assert(boundary_image.shape == laplacian.shape) + est = boundary_image.copy() + + est[1:-1, 1:-1, :] = np.random.random((m-2, n-2, p)) + loss = [] + + for i in range(num_iters): + old_est = est.copy() + est[1:-1, 1:-1, :] = 0.25*(est[0:-2, 1:-1, :] + est[1:-1, 0:-2, :] + + est[2:, 1:-1, :] + est[1:-1, 2:, :] - h*h*laplacian[1:-1, 1:-1, :]) + error = np.sum(np.square(est-old_est)) + loss.append(error) + + return (est) def image_threshold(image, threshold=0.5): - ''' - Threshold the image to make all its elements greater than threshold*MAX = 1 - ''' - m, M = np.min(image), np.max(image) - im = PlotImage(image) - im[im >= threshold] = 1 - im[im < 1] = 0 - return im + ''' + Threshold the image to make all its elements greater than threshold*MAX = 1 + ''' + m, M = np.min(image), np.max(image) + im = PlotImage(image) + im[im >= threshold] = 1 + im[im < 1] = 0 + return im def crop_watermark(gradx, grady, threshold=0.4, boundary_size=2): - """ - Crops the watermark by taking the edge map of magnitude of grad(W) - Assumes the gradx and grady to be in 3 channels - @param: threshold - gives the threshold param - @param: boundary_size - boundary around cropped image - """ - W_mod = np.sqrt(np.square(gradx) + np.square(grady)) - W_mod = PlotImage(W_mod) - W_gray = image_threshold(np.average(W_mod, axis=2), threshold=threshold) - x, y = np.where(W_gray == 1) - - xm, xM = np.min(x) - boundary_size - 1, np.max(x) + boundary_size + 1 - ym, yM = np.min(y) - boundary_size - 1, np.max(y) + boundary_size + 1 - - return gradx[xm:xM, ym:yM, :] , grady[xm:xM, ym:yM, :] + """ + Crops the watermark by taking the edge map of magnitude of grad(W) + Assumes the gradx and grady to be in 3 channels + @param: threshold - gives the threshold param + @param: boundary_size - boundary around cropped image + """ + W_mod = np.sqrt(np.square(gradx) + np.square(grady)) + # Map image values to [0, 1] + W_mod = PlotImage(W_mod) + # Threshold the image with threshold=0.4 + W_gray = image_threshold(np.average(W_mod, axis=2), threshold=threshold) + x, y = np.where(W_gray == 1) + + # Boundary of cropped image (contain watermark) + xm, xM = np.min(x) - boundary_size - 1, np.max(x) + boundary_size + 1 + ym, yM = np.min(y) - boundary_size - 1, np.max(y) + boundary_size + 1 + + return gradx[xm:xM, ym:yM, :], grady[xm:xM, ym:yM, :] def normalized(img): - """ - Return the image between -1 to 1 so that its easier to find out things like - correlation between images, convolutionss, etc. - Currently required for Chamfer distance for template matching. - """ - return (2*PlotImage(img)-1) + """ + Return the image between -1 to 1 so that its easier to find out things like + correlation between images, convolutionss, etc. + Currently required for Chamfer distance for template matching. + """ + return (2*PlotImage(img)-1) + def watermark_detector(img, gx, gy, thresh_low=200, thresh_high=220, printval=False): - """ - Compute a verbose edge map using Canny edge detector, take its magnitude. - Assuming cropped values of gradients are given. - Returns image, start and end coordinates - """ - Wm = (np.average(np.sqrt(np.square(gx) + np.square(gy)), axis=2)) - - img_edgemap = (cv2.Canny(img, thresh_low, thresh_high)) - chamfer_dist = cv2.filter2D(img_edgemap.astype(float), -1, Wm) - - rect = Wm.shape - index = np.unravel_index(np.argmax(chamfer_dist), img.shape[:-1]) - if printval: - print(index) - - x,y = (index[0]-rect[0]/2), (index[1]-rect[1]/2) - im = img.copy() - cv2.rectangle(im, (y, x), (y+rect[1], x+rect[0]), (255, 0, 0)) - return (im, (x, y), (rect[0], rect[1])) + """ + Compute a verbose edge map using Canny edge detector, take its magnitude. + Assuming cropped values of gradients are given. + Returns image, start and end coordinates + """ + Wm = (np.average(np.sqrt(np.square(gx) + np.square(gy)), axis=2)) + + img_edgemap = cv2.Canny(img, thresh_low, thresh_high) + chamfer_dist = cv2.filter2D(img_edgemap.astype(float), -1, Wm) + + rect = Wm.shape + index = np.unravel_index(np.argmax(chamfer_dist), img.shape[:-1]) + if printval: + print(index) + + x = int(index[0]-rect[0]/2) + y = int(index[1]-rect[1]/2) + im = img.copy() + cv2.rectangle(im, (y, x), (y+rect[1], x+rect[0]), (255, 0, 0)) + return (im, (x, y), (rect[0], rect[1])) diff --git a/src/watermark_reconstruct.py b/src/watermark_reconstruct.py index 52aa911..6cce410 100644 --- a/src/watermark_reconstruct.py +++ b/src/watermark_reconstruct.py @@ -4,10 +4,11 @@ import scipy from scipy.sparse import * from scipy.sparse import linalg -from estimate_watermark import * -from closed_form_matting import * +from .estimate_watermark import * +from .closed_form_matting import * from numpy import nan, isnan + def get_cropped_images(foldername, num_images, start, end, shape): ''' This is the part where we get all the images, extract their parts, and then add it to our matrix @@ -32,9 +33,9 @@ def get_cropped_images(foldername, num_images, start, end, shape): _img = _img[_s[0]:(_s[0]+_e[0]), _s[1]:(_s[1]+_e[1]), :] # add to list images images_cropped[index, :, :, :] = _img - index+=1 + index += 1 else: - print("%s not found."%(file)) + print("%s not found." % (file)) return (images_cropped, image_paths) @@ -49,6 +50,8 @@ def _get_ysobel_coord(coord, shape): ] # get sobel coordinates for x + + def _get_xsobel_coord(coord, shape): i, j, k = coord m, n, p = shape @@ -58,14 +61,18 @@ def _get_xsobel_coord(coord, shape): ] # filter + + def _filter_list_item(coord, shape): i, j, k, v = coord m, n, p = shape - if i>=0 and i=0 and j= 0 and i < m and j >= 0 and j < n: return True # Change to ravel index # also filters the wrong guys + + def _change_to_ravel_index(li, shape): li = filter(lambda x: _filter_list_item(x, shape), li) i, j, k, v = zip(*li) @@ -73,13 +80,16 @@ def _change_to_ravel_index(li, shape): # TODO: Consider wrap around of indices to remove the edge at the end of sobel # get Sobel sparse matrix for Y + + def get_ySobel_matrix(m, n, p): size = m*n*p shape = (m, n, p) i, j, k = np.unravel_index(np.arange(size), (m, n, p)) ijk = zip(list(i), list(j), list(k)) ijk_nbrs = map(lambda x: _get_ysobel_coord(x, shape), ijk) - ijk_nbrs_to_index = map(lambda l: _change_to_ravel_index(l, shape), ijk_nbrs) + ijk_nbrs_to_index = map( + lambda l: _change_to_ravel_index(l, shape), ijk_nbrs) # we get a list of idx, values for a particular idx # we have got the complete list now, map it to actual index actual_map = [] @@ -98,7 +108,8 @@ def get_xSobel_matrix(m, n, p): i, j, k = np.unravel_index(np.arange(size), (m, n, p)) ijk = zip(list(i), list(j), list(k)) ijk_nbrs = map(lambda x: _get_xsobel_coord(x, shape), ijk) - ijk_nbrs_to_index = map(lambda l: _change_to_ravel_index(l, shape), ijk_nbrs) + ijk_nbrs_to_index = map( + lambda l: _change_to_ravel_index(l, shape), ijk_nbrs) # we get a list of idx, values for a particular idx # we have got the complete list now, map it to actual index actual_map = [] @@ -110,10 +121,13 @@ def get_xSobel_matrix(m, n, p): return coo_matrix((vals, (i, j)), shape=(size, size)) # get estimated normalized alpha matte -def estimate_normalized_alpha(J, W_m, num_images=30, threshold=170, invert=False, adaptive=False, adaptive_threshold=21, c2=10): + + +def estimate_normalized_alpha(J, W_m, num_images, threshold=170, invert=False, adaptive=False, adaptive_threshold=21, c2=10): _Wm = (255*PlotImage(np.average(W_m, axis=2))).astype(np.uint8) if adaptive: - thr = cv2.adaptiveThreshold(_Wm, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY, adaptive_threshold, c2) + thr = cv2.adaptiveThreshold( + _Wm, 255, cv2.ADAPTIVE_THRESH_GAUSSIAN_C, cv2.THRESH_BINARY, adaptive_threshold, c2) else: ret, thr = cv2.threshold(_Wm, threshold, 255, cv2.THRESH_BINARY) @@ -125,9 +139,9 @@ def estimate_normalized_alpha(J, W_m, num_images=30, threshold=170, invert=False alpha = np.zeros((num_images, m, n)) iterpatch = 900 - print("Estimating normalized alpha using %d images."%(num_images)) + print("Estimating normalized alpha using %d images." % (num_images)) # for all images, calculate alpha - for idx in xrange(num_images): + for idx in range(num_images): imgcopy = thr alph = closed_form_matte(J[idx], imgcopy) alpha[idx] = alph @@ -135,13 +149,14 @@ def estimate_normalized_alpha(J, W_m, num_images=30, threshold=170, invert=False alpha = np.median(alpha, axis=0) return alpha + def estimate_blend_factor(J, W_m, alph, threshold=0.01*255): K, m, n, p = J.shape Jm = (J - W_m) gx_jm = np.zeros(J.shape) gy_jm = np.zeros(J.shape) - for i in xrange(K): + for i in range(K): gx_jm[i] = cv2.Sobel(Jm[i], cv2.CV_64F, 1, 0, 3) gy_jm[i] = cv2.Sobel(Jm[i], cv2.CV_64F, 0, 1, 3) @@ -153,8 +168,9 @@ def estimate_blend_factor(J, W_m, alph, threshold=0.01*255): estIk_grad = np.sqrt(gx_estIk**2 + gy_estIk**2) C = [] - for i in xrange(3): - c_i = np.sum(Jm_grad[:,:,:,i]*estIk_grad[:,:,i])/np.sum(np.square(estIk_grad[:,:,i]))/K + for i in range(3): + c_i = np.sum(Jm_grad[:, :, :, i]*estIk_grad[:, :, i]) / \ + np.sum(np.square(estIk_grad[:, :, i]))/K print(c_i) C.append(c_i) @@ -164,9 +180,11 @@ def estimate_blend_factor(J, W_m, alph, threshold=0.01*255): def Func_Phi(X, epsilon=1e-3): return np.sqrt(X + epsilon**2) + def Func_Phi_deriv(X, epsilon=1e-3): return 0.5/Func_Phi(X, epsilon) + def solve_images(J, W_m, alpha, W_init, gamma=1, beta=1, lambda_w=0.005, lambda_i=1, lambda_a=0.01, iters=4): ''' Master solver, follows the algorithm given in the supplementary. @@ -181,7 +199,7 @@ def solve_images(J, W_m, alpha, W_init, gamma=1, beta=1, lambda_w=0.005, lambda_ sobely = get_ySobel_matrix(m, n, p) Ik = np.zeros(J.shape) Wk = np.zeros(J.shape) - for i in xrange(K): + for i in range(K): Ik[i] = J[i] - W_m Wk[i] = W_init.copy() @@ -189,10 +207,10 @@ def solve_images(J, W_m, alpha, W_init, gamma=1, beta=1, lambda_w=0.005, lambda_ W = W_init.copy() # Iterations - for _ in xrange(iters): + for _ in range(iters): print("------------------------------------") - print("Iteration: %d"%(_)) + print("Iteration: %d" % (_)) # Step 1 print("Step 1") @@ -208,7 +226,7 @@ def solve_images(J, W_m, alpha, W_init, gamma=1, beta=1, lambda_w=0.005, lambda_ alpha_diag = diags(alpha.reshape(-1)) alpha_bar_diag = diags((1-alpha).reshape(-1)) - for i in xrange(K): + for i in range(K): # prep vars Wkx = cv2.Sobel(Wk[i], cv2.CV_64F, 1, 0, 3) Wky = cv2.Sobel(Wk[i], cv2.CV_64F, 0, 1, 3) @@ -218,35 +236,48 @@ def solve_images(J, W_m, alpha, W_init, gamma=1, beta=1, lambda_w=0.005, lambda_ alphaWk = alpha*Wk[i] alphaWk_gx = cv2.Sobel(alphaWk, cv2.CV_64F, 1, 0, 3) - alphaWk_gy = cv2.Sobel(alphaWk, cv2.CV_64F, 0, 1, 3) - - phi_data = diags( Func_Phi_deriv(np.square(alpha*Wk[i] + (1-alpha)*Ik[i] - J[i]).reshape(-1)) ) - phi_W = diags( Func_Phi_deriv(np.square( np.abs(alpha_gx)*Wkx + np.abs(alpha_gy)*Wky ).reshape(-1)) ) - phi_I = diags( Func_Phi_deriv(np.square( np.abs(alpha_gx)*Ikx + np.abs(alpha_gy)*Iky ).reshape(-1)) ) - phi_f = diags( Func_Phi_deriv( ((Wm_gx - alphaWk_gx)**2 + (Wm_gy - alphaWk_gy)**2 ).reshape(-1)) ) - phi_aux = diags( Func_Phi_deriv(np.square(Wk[i] - W).reshape(-1)) ) - phi_rI = diags( Func_Phi_deriv( np.abs(alpha_gx)*(Ikx**2) + np.abs(alpha_gy)*(Iky**2) ).reshape(-1) ) - phi_rW = diags( Func_Phi_deriv( np.abs(alpha_gx)*(Wkx**2) + np.abs(alpha_gy)*(Wky**2) ).reshape(-1) ) - - L_i = sobelx.T.dot(cx*phi_rI).dot(sobelx) + sobely.T.dot(cy*phi_rI).dot(sobely) - L_w = sobelx.T.dot(cx*phi_rW).dot(sobelx) + sobely.T.dot(cy*phi_rW).dot(sobely) - L_f = sobelx.T.dot(phi_f).dot(sobelx) + sobely.T.dot(phi_f).dot(sobely) + alphaWk_gy = cv2.Sobel(alphaWk, cv2.CV_64F, 0, 1, 3) + + phi_data = diags(Func_Phi_deriv( + np.square(alpha*Wk[i] + (1-alpha)*Ik[i] - J[i]).reshape(-1))) + phi_W = diags(Func_Phi_deriv( + np.square(np.abs(alpha_gx)*Wkx + np.abs(alpha_gy)*Wky).reshape(-1))) + phi_I = diags(Func_Phi_deriv( + np.square(np.abs(alpha_gx)*Ikx + np.abs(alpha_gy)*Iky).reshape(-1))) + phi_f = diags(Func_Phi_deriv( + ((Wm_gx - alphaWk_gx)**2 + (Wm_gy - alphaWk_gy)**2).reshape(-1))) + phi_aux = diags(Func_Phi_deriv(np.square(Wk[i] - W).reshape(-1))) + phi_rI = diags(Func_Phi_deriv(np.abs(alpha_gx) * + (Ikx**2) + np.abs(alpha_gy)*(Iky**2)).reshape(-1)) + phi_rW = diags(Func_Phi_deriv(np.abs(alpha_gx) * + (Wkx**2) + np.abs(alpha_gy)*(Wky**2)).reshape(-1)) + + L_i = sobelx.T.dot(cx*phi_rI).dot(sobelx) + \ + sobely.T.dot(cy*phi_rI).dot(sobely) + L_w = sobelx.T.dot(cx*phi_rW).dot(sobelx) + \ + sobely.T.dot(cy*phi_rW).dot(sobely) + L_f = sobelx.T.dot(phi_f).dot(sobelx) + \ + sobely.T.dot(phi_f).dot(sobely) A_f = alpha_diag.T.dot(L_f).dot(alpha_diag) + gamma*phi_aux - bW = alpha_diag.dot(phi_data).dot(J[i].reshape(-1)) + beta*L_f.dot(W_m.reshape(-1)) + gamma*phi_aux.dot(W.reshape(-1)) + bW = alpha_diag.dot(phi_data).dot( + J[i].reshape(-1)) + beta*L_f.dot(W_m.reshape(-1)) + gamma*phi_aux.dot(W.reshape(-1)) bI = alpha_bar_diag.dot(phi_data).dot(J[i].reshape(-1)) - A = vstack([hstack([(alpha_diag**2)*phi_data + lambda_w*L_w + beta*A_f, alpha_diag*alpha_bar_diag*phi_data]), \ - hstack([alpha_diag*alpha_bar_diag*phi_data, (alpha_bar_diag**2)*phi_data + lambda_i*L_i])]).tocsr() + A = vstack([hstack([(alpha_diag**2)*phi_data + lambda_w*L_w + beta*A_f, alpha_diag*alpha_bar_diag*phi_data]), + hstack([alpha_diag*alpha_bar_diag*phi_data, (alpha_bar_diag**2)*phi_data + lambda_i*L_i])]).tocsr() b = np.hstack([bW, bI]) x = linalg.spsolve(A, b) - + Wk[i] = x[:size].reshape(m, n, p) Ik[i] = x[size:].reshape(m, n, p) - plt.subplot(3,1,1); plt.imshow(PlotImage(J[i])) - plt.subplot(3,1,2); plt.imshow(PlotImage(Wk[i])) - plt.subplot(3,1,3); plt.imshow(PlotImage(Ik[i])) + plt.subplot(3, 1, 1) + plt.imshow(PlotImage(J[i])) + plt.subplot(3, 1, 2) + plt.imshow(PlotImage(Wk[i])) + plt.subplot(3, 1, 3) + plt.imshow(PlotImage(Ik[i])) plt.draw() plt.pause(0.001) print(i) @@ -258,39 +289,45 @@ def solve_images(J, W_m, alpha, W_init, gamma=1, beta=1, lambda_w=0.005, lambda_ plt.imshow(PlotImage(W)) plt.draw() plt.pause(0.001) - + # Step 3 print("Step 3") W_diag = diags(W.reshape(-1)) - + for i in range(K): alphaWk = alpha*Wk[i] alphaWk_gx = cv2.Sobel(alphaWk, cv2.CV_64F, 1, 0, 3) - alphaWk_gy = cv2.Sobel(alphaWk, cv2.CV_64F, 0, 1, 3) - phi_f = diags( Func_Phi_deriv( ((Wm_gx - alphaWk_gx)**2 + (Wm_gy - alphaWk_gy)**2 ).reshape(-1)) ) - - phi_kA = diags(( (Func_Phi_deriv((((alpha*Wk[i] + (1-alpha)*Ik[i] - J[i])**2)))) * ((W-Ik[i])**2) ).reshape(-1)) - phi_kB = (( (Func_Phi_deriv((((alpha*Wk[i] + (1-alpha)*Ik[i] - J[i])**2))))*(W-Ik[i])*(J[i]-Ik[i]) ).reshape(-1)) - - phi_alpha = diags(Func_Phi_deriv(alpha_gx**2 + alpha_gy**2).reshape(-1)) - L_alpha = sobelx.T.dot(phi_alpha.dot(sobelx)) + sobely.T.dot(phi_alpha.dot(sobely)) - - L_f = sobelx.T.dot(phi_f).dot(sobelx) + sobely.T.dot(phi_f).dot(sobely) + alphaWk_gy = cv2.Sobel(alphaWk, cv2.CV_64F, 0, 1, 3) + phi_f = diags(Func_Phi_deriv( + ((Wm_gx - alphaWk_gx)**2 + (Wm_gy - alphaWk_gy)**2).reshape(-1))) + + phi_kA = diags(((Func_Phi_deriv( + (((alpha*Wk[i] + (1-alpha)*Ik[i] - J[i])**2)))) * ((W-Ik[i])**2)).reshape(-1)) + phi_kB = (((Func_Phi_deriv( + (((alpha*Wk[i] + (1-alpha)*Ik[i] - J[i])**2))))*(W-Ik[i])*(J[i]-Ik[i])).reshape(-1)) + + phi_alpha = diags(Func_Phi_deriv( + alpha_gx**2 + alpha_gy**2).reshape(-1)) + L_alpha = sobelx.T.dot(phi_alpha.dot(sobelx)) + \ + sobely.T.dot(phi_alpha.dot(sobely)) + + L_f = sobelx.T.dot(phi_f).dot(sobelx) + \ + sobely.T.dot(phi_f).dot(sobely) A_tilde_f = W_diag.T.dot(L_f).dot(W_diag) # Ax = b, setting up A - if i==0: + if i == 0: A1 = phi_kA + lambda_a*L_alpha + beta*A_tilde_f b1 = phi_kB + beta*W_diag.dot(L_f).dot(W_m.reshape(-1)) else: A1 += (phi_kA + lambda_a*L_alpha + beta*A_tilde_f) b1 += (phi_kB + beta*W_diag.T.dot(L_f).dot(W_m.reshape(-1))) - alpha = linalg.spsolve(A1, b1).reshape(m,n,p) + alpha = linalg.spsolve(A1, b1).reshape(m, n, p) plt.imshow(PlotImage(alpha)) plt.draw() plt.pause(0.001) - + return (Wk, Ik, W, alpha) @@ -299,7 +336,7 @@ def changeContrastImage(J, I): cJ2 = J[-1, -1, :] cI1 = I[0, 0, :] - cI2 = I[-1,-1, :] + cI2 = I[-1, -1, :] I_m = cJ1 + (I-cI1)/(cI2-cI1)*(cJ2-cJ1) - return I_m \ No newline at end of file + return I_m