| 7 | 8 |
new file mode 100755 |
| ... | ... |
@@ -0,0 +1,89 @@ |
| 1 |
+#!/usr/bin/env python3 |
|
| 2 |
+ |
|
| 3 |
+import sys |
|
| 4 |
+import PIL |
|
| 5 |
+import numpy as np |
|
| 6 |
+import scipy.signal |
|
| 7 |
+import scipy.ndimage |
|
| 8 |
+import matplotlib.colors |
|
| 9 |
+import matplotlib.pyplot as plt |
|
| 10 |
+ |
|
| 11 |
+INPUT = sys.argv[1] |
|
| 12 |
+OUTPUT = len(sys.argv) > 2 and sys.argv[2] or None |
|
| 13 |
+HIGHLIGHT = [255, 0, 0] |
|
| 14 |
+ |
|
| 15 |
+subplot_index = 0 |
|
| 16 |
+def show(title, *args, **kwargs): |
|
| 17 |
+ global subplot_index |
|
| 18 |
+ subplot_index += 1 |
|
| 19 |
+ plt.subplot(2, 3, subplot_index) |
|
| 20 |
+ plt.imshow(*args, **kwargs) |
|
| 21 |
+ plt.title(title) |
|
| 22 |
+ |
|
| 23 |
+def blur(im, sigma, cutoff=0.01): |
|
| 24 |
+ length = int(2 * sigma * np.sqrt(-2*np.log(cutoff))) |
|
| 25 |
+ signal = scipy.signal.windows.gaussian(length, sigma) |
|
| 26 |
+ kernel = np.outer(signal, signal) / sum(signal)**2 |
|
| 27 |
+ return scipy.signal.fftconvolve(im, kernel, mode='same') |
|
| 28 |
+ |
|
| 29 |
+def vec(x, y, adjust=False): |
|
| 30 |
+ mag = np.sqrt(x*x + y*y) |
|
| 31 |
+ mag /= np.max(mag) |
|
| 32 |
+ ang = np.arctan2(y, x) / (2*np.pi) + 0.5 |
|
| 33 |
+ if adjust: |
|
| 34 |
+ min = np.min(ang) |
|
| 35 |
+ max = np.max(ang) |
|
| 36 |
+ ang = (ang - min) / (max - min) |
|
| 37 |
+ return matplotlib.colors.hsv_to_rgb(np.moveaxis((ang, mag, mag), 0, 2)) |
|
| 38 |
+ |
|
| 39 |
+# Load image |
|
| 40 |
+gray = np.array(PIL.Image.open(INPUT).convert('L'), dtype='float32')
|
|
| 41 |
+sigma = np.hypot(*gray.shape) / 30 |
|
| 42 |
+ |
|
| 43 |
+# Calculate gradients |
|
| 44 |
+dy, dx = np.gradient(gray) |
|
| 45 |
+show("Gradient", vec(dx, dy))
|
|
| 46 |
+ |
|
| 47 |
+# Fold gradients |
|
| 48 |
+sel = (dx if np.sum(np.abs(dx)) > np.sum(np.abs(dy)) else dy) < 0 |
|
| 49 |
+dx[sel] = -dx[sel] |
|
| 50 |
+dy[sel] = -dy[sel] |
|
| 51 |
+show("Gradient, folded", vec(dx, dy))
|
|
| 52 |
+ |
|
| 53 |
+# Blur gradients |
|
| 54 |
+dx = blur(dx, sigma) |
|
| 55 |
+dy = blur(dy, sigma) |
|
| 56 |
+show("Gradient, blurred", vec(dx, dy, True))
|
|
| 57 |
+ |
|
| 58 |
+# Calculate dot product with maximal blurred gradient, normalized |
|
| 59 |
+mag = dx*dx + dy*dy |
|
| 60 |
+ind = np.unravel_index(np.argmax(mag), mag.shape) |
|
| 61 |
+mdx = dx[ind] |
|
| 62 |
+mdy = dy[ind] |
|
| 63 |
+weight = (dx*mdx + dy*mdy) / (mdx*mdx + mdy*mdy) |
|
| 64 |
+show("Weight", weight)
|
|
| 65 |
+ |
|
| 66 |
+# Mask and display |
|
| 67 |
+mask = weight > 2/3 |
|
| 68 |
+im = np.array(PIL.Image.open(INPUT)) |
|
| 69 |
+im[mask] = (im[mask] + HIGHLIGHT) / 2 |
|
| 70 |
+show("Mask", im)
|
|
| 71 |
+ |
|
| 72 |
+# Rotate and crop |
|
| 73 |
+angle = np.degrees(np.arctan2(mdy, mdx)) |
|
| 74 |
+gray = scipy.ndimage.rotate(gray, angle) |
|
| 75 |
+mask = scipy.ndimage.rotate(mask, angle) |
|
| 76 |
+my, mx = np.nonzero(mask) |
|
| 77 |
+mx = np.min(mx), np.max(mx) |
|
| 78 |
+my = np.min(my), np.max(my) |
|
| 79 |
+cx = int((mx[1] - mx[0]) * 0.2) |
|
| 80 |
+cy = int((my[1] - my[0]) * 0.3) |
|
| 81 |
+crop = gray[my[0]+cy:my[1]-cy, mx[0]-cx:mx[1]+cx] |
|
| 82 |
+show("Rotated, cropped", crop, 'gray')
|
|
| 83 |
+ |
|
| 84 |
+# Save/show |
|
| 85 |
+plt.tight_layout() |
|
| 86 |
+if OUTPUT: |
|
| 87 |
+ plt.savefig(OUTPUT, dpi=200) |
|
| 88 |
+else: |
|
| 89 |
+ plt.show() |
git.rcrnstn.net
