import imageio
from skimage.color import rgb2gray
from skimage import transform
img = imageio.imread("../raw/ZIMBABWE/EVAP.jpeg")
/tmp/ipykernel_355275/3550911579.py:1: DeprecationWarning: Starting with ImageIO v3 the behavior of this function will switch to that of iio.v3.imread. To keep the current behavior (and make this warning disappear) use `import imageio.v2 as imageio` or call `imageio.v2.imread` directly.
img = imageio.imread("../raw/ZIMBABWE/EVAP.jpeg")
import matplotlib.pyplot as plt
plt.imshow(img)
<matplotlib.image.AxesImage at 0x7f25ee0f9520>
img_gray = rgb2gray(img)
import seaborn as sns
sns.histplot(img_gray.flatten())
<Axes: ylabel='Count'>
img_gray.std()
0.2228197984558268
from skimage.feature import canny
from random import choice
import numpy as np
def detect_lines(img_gray, plot=False):
edges = canny(img_gray, sigma=img_gray.std(), low_threshold=0.5, high_threshold=0.98)
lines = transform.probabilistic_hough_line(edges, line_length=150, line_gap=10)
if plot:
fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(ncols=2, nrows=2, figsize=(15,12), sharey=True, dpi=150)
ax0.imshow(img_gray, cmap="gray")
ax1.imshow(edges, cmap="gray")
source_vectors = []
dest_vectors = []
for line in lines:
vector = np.array(line)
xs = vector[:,0]
ys = vector[:,1]
slope = (ys[0] - ys[1]) / (xs[0] - xs[1])
vertical, horizontal, orthogonal = compute_orientation(slope)
if not(vertical is horizontal) and not orthogonal:
source_vectors.append(vector)
vector_straight = np.copy(vector)
if horizontal:
vector_straight[:,1]= vector[:,1].mean()
else:
vector_straight[:,0]= vector[:,0].mean()
dest_vectors.append(vector_straight)
else:
pass
#print("Discard", vector)
if plot:
ax2.imshow(edges*0)
for line in source_vectors:
p0, p1= line
ax2.plot((p0[0], p1[0]), (p0[1], p1[1]))
ax3.imshow(edges*0)
for line in dest_vectors:
p0, p1 = line
ax3.plot((p0[0], p1[0]), (p0[1], p1[1]))
return source_vectors, dest_vectors
def compute_orientation(slope, vtol=0.1, htol=0.05):
angle = np.abs(np.arctan(slope))
vertical = np.isclose(angle, np.pi/2, atol=vtol)
horizontal = np.isclose(angle, 0, atol=htol)
orthogonal = (angle == 0) or (angle == np.pi/2)
return vertical,horizontal,orthogonal
source_vectors, dest_vectors = detect_lines(img_gray, plot=True)
/tmp/ipykernel_355275/1315240660.py:24: RuntimeWarning: divide by zero encountered in scalar divide
slope = (ys[0] - ys[1]) / (xs[0] - xs[1])
import cv2 as cv
from skimage.exposure import rescale_intensity
def detect_contours(img_gray, plot=False):
edges = cv.Canny(rescale_intensity(img_gray, out_range=(0, 255)).astype('uint8'), 128, 250)
# edges = canny(img_gray, sigma=img_gray.std(), low_threshold=0.5, high_threshold=0.98)
contours, _ = cv.findContours(edges, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE)
if plot:
fig, ((ax0, ax1), (ax2, ax3)) = plt.subplots(ncols=2, nrows=2, figsize=(15,12), sharey=True, dpi=150)
ax0.imshow(img_gray, cmap="gray")
ax1.imshow(edges, cmap="gray")
ax2.imshow(edges*0)
ax3.imshow(edges*0)
source_vectors = []
dest_vectors = []
for contour in contours:
xs = contour[...,0]
ys = contour[..., 1]
x_span = xs.max() - xs.min()
y_span = ys.max() - ys.min()
if max(x_span, y_span) < 500:
continue
slope = y_span / x_span
vertical, horizontal, orthogonal = compute_orientation(slope, htol=0.1)
vector = contour
vector_straight = np.copy(contour)
if horizontal:
vector_straight[...,1]= np.median(vector[...,1])
elif vertical:
vector_straight[...,0]= np.median(vector[...,0])
else:
continue
source_vectors.append(vector[:,0,:])
dest_vectors.append(vector_straight[:,0,:])
if plot:
ax2.plot(xs, ys)
ax3.plot(vector_straight[...,0], vector_straight[..., 1])
return source_vectors, dest_vectors
source_vectors, dest_vectors = detect_contours(img_gray, plot=True)
deskew = transform.ProjectiveTransform()
# deskew = transform.PiecewiseAffineTransform()
dest, source = map(np.concatenate, (dest_vectors, source_vectors))
nb_pts = None
deskew.estimate(dest[:nb_pts], source[:nb_pts])
warped = transform.warp(img_gray, deskew)
def compare_warp(im1, im1_anchors, im2, im2_anchors):
fig, ax = plt.subplots(ncols=2, figsize=(10, 5), dpi=200)
ax[0].imshow(im1, cmap=plt.cm.gray)
ax[0].plot(im1_anchors[:nb_pts, 0], im1_anchors[:nb_pts, 1], '.r')
ax[1].imshow(im2, cmap=plt.cm.gray)
ax[1].plot(im2_anchors[:nb_pts, 0], im2_anchors[:nb_pts, 1], '.b')
for a in ax:
a.axis('off')
plt.tight_layout()
plt.show()
return ax
ax = compare_warp(img_gray, source, warped, dest)
dest.shape
(4816, 2)
import cv2 as cv
from random import sample
rows,cols = img_gray.shape
def opencv_compat(arr):
return np.ascontiguousarray(arr, dtype='float32')
# from itertools import islice
# def batched(iterable, n=4):
# # batched('ABCDEFG', 3) --> ABC DEF G
# if n < 1:
# raise ValueError('n must be at least one')
# it = iter(iterable)
# while batch := tuple(islice(it, n)):
# yield batch
# warped = img_gray
# for n, (dest_batch, source_batch) in enumerate(zip(*map(batched, (dest, source)))):
# dest_mat, source_mat = map(opencv_compat, (dest_batch, source_batch))
# dest_mat, source_mat = map(opencv_compat, (dest_batch, source_batch))
# M = cv.getAffineTransform(dest_mat[:3], source_mat[:3])
# warped = cv.warpAffine(warped, M, (cols,rows))
# if n == 2:
# break
dest_mat, source_mat = map(opencv_compat, (dest, source))
subset_mat = sample(range(dest_mat.shape[0]), 4)
M = cv.getPerspectiveTransform(dest_mat[subset_mat], source_mat[subset_mat])
print(M)
warped = cv.warpPerspective(img_gray, M, (cols,rows))
M, mask = cv.findHomography(dest, source, method=cv.RANSAC)
# subset_mat = sample(range(dest_mat.shape[0]), 3)
# M = cv.getAffineTransform(dest_mat[subset_mat], source_mat[subset_mat])
# print(M)
# warped = cv.warpAffine(img_gray, M, (cols,rows))
_ = compare_warp(img_gray, source_mat[subset_mat], warped, dest_mat[subset_mat])
[[ 1.03600239e+00 -2.37998114e-03 -2.08099663e+01]
[-4.72896490e-03 9.88762130e-01 4.17156934e+00]
[-7.76485516e-06 -4.50478530e-06 1.00000000e+00]]
dest.shape, source.shape
((4816, 2), (4816, 2))
import scipy
from tps import ThinPlateSpline
# Create the tps object
tps = ThinPlateSpline(alpha=0.5) # Regularization
# Fit the control and target points
tps.fit(source, dest)
width, height = img_gray.shape
# Create the 2d meshgrid of indices for output image
output_indices = np.indices((height, width), dtype=np.float64).transpose(1, 2, 0) # Shape: (H, W, 2)
# Transform it into the input indices
input_indices = tps.transform(output_indices.reshape(-1, 2)).reshape(output_indices.shape)
# Interpolate the resulting image
warped = scipy.ndimage.map_coordinates(img_gray, input_indices.transpose(2, 0, 1))
ax = compare_warp(img_gray, source, warped, np.zeros([2,2]))
---------------------------------------------------------------------------
MemoryError Traceback (most recent call last)
/local_disk/data/ai-for-obs/notebooks/2023-11-03-deskew.ipynb Cell 13 line 1
<a href='vscode-notebook-cell:/local_disk/data/ai-for-obs/notebooks/2023-11-03-deskew.ipynb#X14sZmlsZQ%3D%3D?line=12'>13</a> output_indices = np.indices((height, width), dtype=np.float64).transpose(1, 2, 0) # Shape: (H, W, 2)
<a href='vscode-notebook-cell:/local_disk/data/ai-for-obs/notebooks/2023-11-03-deskew.ipynb#X14sZmlsZQ%3D%3D?line=14'>15</a> # Transform it into the input indices
---> <a href='vscode-notebook-cell:/local_disk/data/ai-for-obs/notebooks/2023-11-03-deskew.ipynb#X14sZmlsZQ%3D%3D?line=15'>16</a> input_indices = tps.transform(output_indices.reshape(-1, 2)).reshape(output_indices.shape)
<a href='vscode-notebook-cell:/local_disk/data/ai-for-obs/notebooks/2023-11-03-deskew.ipynb#X14sZmlsZQ%3D%3D?line=17'>18</a> # Interpolate the resulting image
<a href='vscode-notebook-cell:/local_disk/data/ai-for-obs/notebooks/2023-11-03-deskew.ipynb#X14sZmlsZQ%3D%3D?line=18'>19</a> warped = scipy.ndimage.map_coordinates(img_gray, input_indices.transpose(2, 0, 1))
File ~/.pyenv/versions/3.9.13/envs/ai-for-obs/lib/python3.9/site-packages/tps/thin_plate_spline.py:110, in ThinPlateSpline.transform(self, X)
107 X = _ensure_2d(X)
108 assert X.shape[1] == self.control_points.shape[1]
--> 110 phi = self._radial_distance(X) # n x n_c
112 X = np.hstack([phi, np.ones((X.shape[0], 1)), X]) # n x (n_c + 1 + d_s)
113 return X @ self.parameters
File ~/.pyenv/versions/3.9.13/envs/ai-for-obs/lib/python3.9/site-packages/tps/thin_plate_spline.py:128, in ThinPlateSpline._radial_distance(self, X)
115 def _radial_distance(self, X: np.ndarray) -> np.ndarray:
116 """Compute the pairwise radial distances of the given points to the control points
117
118 Input dimensions are not checked.
(...)
126 Shape: (n, n_c)
127 """
--> 128 dist = cdist(X, self.control_points)
129 dist[dist == 0] = 1 # phi(r) = r^2 log(r) -> (phi(0) = 0)
130 return dist**2 * np.log(dist)
File ~/.pyenv/versions/3.9.13/envs/ai-for-obs/lib/python3.9/site-packages/scipy/spatial/distance.py:3006, in cdist(XA, XB, metric, out, **kwargs)
3004 cdist_fn = metric_info.cdist_func
3005 _extra_windows_error_checks(XA, out, (mA, mB), **kwargs)
-> 3006 return cdist_fn(XA, XB, out=out, **kwargs)
3007 elif mstr.startswith("test_"):
3008 metric_info = _TEST_METRICS.get(mstr, None)
MemoryError: Unable to allocate 139. GiB for an array with shape (3869030, 4816) and data type float64
output_indices.shape
(1870, 2069, 2)