# Copyright (c) 2016-2017 Weitian LI # MIT license """ Image (only gray-scale image, i.e., matrix) transformation utilities. References ---------- - Leptonica: Rotation http://www.leptonica.com/rotation.html - Image rotation by MATLAB without using imrotate https://stackoverflow.com/a/19687481/4856091 https://stackoverflow.com/a/19689081/4856091 - Stackoverflow: Python: Rotating greyscale images https://codereview.stackexchange.com/a/41903 """ import numpy as np import numba as nb from scipy import ndimage @nb.jit([nb.float64[:, :](nb.int64[:, :], nb.float64, nb.boolean, nb.boolean, nb.float64), nb.float64[:, :](nb.float64[:, :], nb.float64, nb.boolean, nb.boolean, nb.float64)], nopython=True) def rotate_center(imgin, angle, interp=True, reshape=True, fill_value=0.0): """ Rotate the input image (only gray-scale image currently supported) by a given angle about its center. Parameters ---------- imgin : 2D `~numpy.ndarray` Input gray-scale image to be rotated angle : float Rotation angle (unit: [ degree ]) interp : bool, optional Use the area mapping of the 4 closest input pixels (``interp=True``), which is also the same as "bilinear interpolation", or use the nearest neighbor (``interp=False``) for rotated pixels. reshape : bool, optional Whether adapt the output shape so that the input image is contained completely in the output? fill_value : float, optional Value used to fill pixels in the rotated image that corresponding pixels outside the boundaries of the input image. """ nrow, ncol = imgin.shape # Rotation transformation image angle = np.deg2rad(angle) mrotate = np.zeros((2, 2), dtype=np.float64) mrotate[0, 0] = np.cos(angle) mrotate[0, 1] = np.sin(angle) mrotate[1, 0] = -np.sin(angle) mrotate[1, 1] = np.cos(angle) # Determine the shape of rotated image corner00 = np.array((0, 0)) corner01 = np.array((0, ncol-1)) corner10 = np.array((nrow-1, 0)) corner11 = np.array((nrow-1, ncol-1)) corners = np.vstack((corner00, corner01, corner10, corner11)) if reshape: dest = np.dot(corners.astype(np.float64), mrotate) # XXX: ``numba`` does not support ``np.max()`` with arguments minr = np.min(dest[:, 0]) minc = np.min(dest[:, 1]) maxr = np.max(dest[:, 0]) maxc = np.max(dest[:, 1]) nr = int(maxr - minr + 0.5) nc = int(maxc - minc + 0.5) else: # Constraint to be same shape nr = nrow nc = ncol imgout = np.ones((nr, nc)) * fill_value # # Calculate the offset, for easier transformation of rotated pixels # back to input image. # # NOTE: # Notations: # P_out : (r_out, c_out) a pixel in the output rotated image # C_out : center position of the output rotated image # P_in : (r_in, c_in) a pixel in the input image # C_in : center position of the input image # R : rotation matrix # R_T : transposed rotation matrix # The rotation relation between pixel pair is (about the center): # (P_in - C_in) * R = P_out - C_out # Then: # (P_in - C_in) = (P_out - C_out) * R_T # And then: # P_in = C_in + (P_out-C_out) * R_T = P_out*R_T + (C_in - C_out*R_T) # Thus can define the "offset" as: # offset = C_in - C_out * R_T # Then the transformation back to input image is simply given by: # P_in = P_out * R_T + offset # center_in = np.array((nrow/2.0 - 0.5, ncol/2.0 - 0.5)) center_out = np.array((nr/2.0 - 0.5, nc/2.0 - 0.5)) mrotate_T = mrotate.transpose() offset = center_in - np.dot(center_out, mrotate_T) # Map pixels of the rotated image to the input image for rr in range(nr): for cc in range(nc): p_out = np.array((rr, cc)) p_in = np.dot(p_out.astype(np.float64), mrotate_T) + offset if np.all((p_in > corner00) & (p_in < corner11)): # Calculate the pixel value for the rotated image if interp: # Use area mapping of the 4 closest input pixels idx_rf, idx_cf = np.floor(p_in).astype(np.int64) idx_rc, idx_cc = np.ceil(p_in).astype(np.int64) # NOTE: # It is possible that ``p_in[0]`` and/or ``p_in[1]`` # are just integers, which cause ``idx_rf == idx_rc`` # and/or ``idx_cf == idx_cc``, which further lead to # the calculated pixel value ``p_val = 0``. if idx_rf == idx_rc: idx_rc += 1 if idx_cf == idx_cc: idx_cc += 1 # Calculate the overlapping areas p_r, p_c = p_in p4_area = np.array([(idx_rc - p_r) * (idx_cc - p_c), (idx_rc - p_r) * (p_c - idx_cf), (p_r - idx_rf) * (idx_cc - p_c), (p_r - idx_rf) * (p_c - idx_cf)]) p4_val = np.array((imgin[idx_rf, idx_cf], imgin[idx_rf, idx_cc], imgin[idx_rc, idx_cf], imgin[idx_rc, idx_cc])) p_val = np.sum(p4_area * p4_val) else: # Use the nearest neighbor as the rotated value idx_r = round(p_in[0]) idx_c = round(p_in[1]) p_val = imgin[idx_r, idx_c] # imgout[rr, cc] = p_val return imgout def circle2ellipse(imgcirc, bfraction, rotation=0.0): """ Shrink the input circle image with respect to the center along the column (axis) to transform the circle to an ellipse, and then rotate around the image center. Parameters ---------- imgcirc : 2D `~numpy.ndarray` Input image grid containing a circle at the center bfraction : float The fraction of the semi-minor axis w.r.t. the semi-major axis (assumed to be the width of the input image), to determine the shrunk size (height) of the output image. Should be a fraction within [0, 1] rotation : float, optional Rotation angle (unit: [ degree ]) Returns ------- imgout : 2D `~numpy.ndarray` Image of the same size as the input circle image. """ nrow, ncol = imgcirc.shape # Shrink the circle to be elliptical nrow2 = nrow * bfraction nrow2 = int(nrow2 / 2) * 2 + 1 # be odd img2 = ndimage.zoom(imgcirc, zoom=(nrow2/nrow, 1.0), order=1) # Pad the shrunk image to have the same size as input imgout = np.zeros(shape=(nrow, ncol)) r1 = int((nrow - nrow2) / 2) r2 = r1 + nrow2 imgout[r1:r2, :] = img2 # Rotate the ellipse imgout = ndimage.rotate(imgout, angle=rotation, reshape=False, order=1) return imgout