Share bilinear interpolation function for other code

Make bilinear_interpolation function callable by other cython code and change
fast_homography to use this function. Also, improve performance of some
functions by making them inline and fix broken test of fast_homography.
This commit is contained in:
Johannes Schönberger
2012-08-16 19:08:13 +02:00
parent 75c7926412
commit 0ca7933a7d
3 changed files with 72 additions and 75 deletions
+12
View File
@@ -0,0 +1,12 @@
cimport numpy as np
import numpy as np
cdef inline double bilinear_interpolation(double* image, int rows, int cols,
double r, double c, char mode,
double cval=*)
cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
char mode, double cval=*)
cdef inline int coord_map(int dim, int coord, char mode)
+49 -63
View File
@@ -1,38 +1,50 @@
#cython: cdivison=True boundscheck=False
#cython: cdivison=True
#cython: boundscheck=False
#cython: nonecheck=False
#cython: wraparound=False
__all__ = ['homography']
cimport cython
cimport numpy as np
import numpy as np
import cython
from cython.operator import dereference
from libc.math cimport ceil, floor
np.import_array()
cdef extern from "math.h":
double floor(double)
double fmod(double, double)
cdef inline double bilinear_interpolation(double* image, int rows, int cols,
double r, double c, char mode,
double cval=0):
cdef double dr, dc
cdef int minr, minc, maxr, maxc
cdef double get_pixel(double *image, int rows, int cols,
int r, int c, char mode, double cval=0):
minr = <int>floor(r)
minc = <int>floor(c)
maxr = <int>ceil(r)
maxc = <int>ceil(c)
dr = r - minr
dc = c - minc
top = (1 - dc) * get_pixel(image, rows, cols, minr, minc, mode, cval) \
+ dc * get_pixel(image, rows, cols, minr, maxc, mode, cval)
bottom = (1 - dc) * get_pixel(image, rows, cols, maxr, minc, mode, cval) \
+ dc * get_pixel(image, rows, cols, maxr, maxc, mode, cval)
return (1 - dr) * top + dr * bottom
cdef inline double get_pixel(double* image, int rows, int cols, int r, int c,
char mode, double cval=0):
"""Get a pixel from the image, taking wrapping mode into consideration.
Parameters
----------
image : *double
image : array of dtype double
Input image.
rows, cols : int
Dimensions of image.
rows, cols: int
Shape of image.
r, c : int
Position at which to get the pixel.
mode : {'C', 'W', 'M'}
Wrapping mode. Constant, Wrap or Mirror.
Wrapping mode. Constant, Wrap or Mirror.
cval : double
Constant value to use for mode constant.
Constant value to use for constant mode.
"""
if mode == 'C':
if (r < 0) or (r > rows - 1) or (c < 0) or (c > cols - 1):
@@ -40,13 +52,13 @@ cdef double get_pixel(double *image, int rows, int cols,
else:
return image[r * cols + c]
else:
return image[coord_map(rows, r, mode) * cols +
coord_map(cols, c, mode)]
return image[coord_map(rows, r, mode) * cols + coord_map(cols, c, mode)]
cdef int coord_map(int dim, int coord, char mode):
cdef inline int coord_map(int dim, int coord, char mode):
"""
Wrap a coordinate, according to a given dimension and mode.
Wrap a coordinate, according to a given mode.
Parameters
----------
dim : int
@@ -56,7 +68,7 @@ cdef int coord_map(int dim, int coord, char mode):
mode : {'W', 'M'}
Whether to wrap or mirror the coordinate if it
falls outside [0, dim).
"""
dim = dim - 1
if mode == 'M': # mirror
@@ -79,7 +91,8 @@ cdef int coord_map(int dim, int coord, char mode):
return coord
cdef tf(double x, double y, double* H, double *x_, double *y_):
cdef inline tf(double x, double y, double* H, double *x_, double *y_):
"""Apply a homography to a coordinate.
Parameters
@@ -98,18 +111,15 @@ cdef tf(double x, double y, double* H, double *x_, double *y_):
yy = H[3] * x + H[4] * y + H[5]
zz = H[6] * x + H[7] * y + H[8]
xx = xx / zz
yy = yy / zz
x_[0] = xx / zz
y_[0] = yy / zz
x_[0] = xx
y_[0] = yy
@cython.boundscheck(False)
def homography(np.ndarray image, np.ndarray H, output_shape=None,
mode='constant', double cval=0):
"""
Projective transformation (homography).
Perform a projective transformation (homography) of a
floating point image, using bi-linear interpolation.
@@ -140,8 +150,6 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
Transformation matrix H that defines the homography.
output_shape : tuple (rows, cols)
Shape of the output image generated.
order : int
Order of splines used in interpolation.
mode : {'constant', 'mirror', 'wrap'}
How to handle values outside the image borders.
cval : string
@@ -150,8 +158,7 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
"""
cdef np.ndarray[dtype=np.double_t, ndim=2, mode="c"] img = \
np.ascontiguousarray(image, dtype=np.double)
cdef np.ndarray[dtype=np.double_t, ndim=2] img = image.astype(np.double)
cdef np.ndarray[dtype=np.double_t, ndim=2, mode="c"] M = \
np.ascontiguousarray(np.linalg.inv(H))
@@ -165,7 +172,6 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
elif mode == 'mirror':
mode_c = ord('M')
cdef int out_r, out_c, columns, rows
if output_shape is None:
out_r = img.shape[0]
out_c = img.shape[1]
@@ -173,37 +179,17 @@ def homography(np.ndarray image, np.ndarray H, output_shape=None,
out_r = output_shape[0]
out_c = output_shape[1]
rows = img.shape[0]
columns = img.shape[1]
cdef np.ndarray[dtype=np.double_t, ndim=2] out = \
np.zeros((out_r, out_c), dtype=np.double)
cdef int tfr, tfc, r_int, c_int
cdef double y0, y1, y2, y3
cdef double r, c, z, t, u
cdef int tfr, tfc
cdef double r, c
cdef int rows = img.shape[0]
cdef int cols = img.shape[1]
for tfr in range(out_r):
for tfc in range(out_c):
tf(tfc, tfr, <double*>M.data, &c, &r)
r_int = <int>floor(r)
c_int = <int>floor(c)
t = r - r_int
u = c - c_int
y0 = get_pixel(<double*>img.data, rows, columns,
r_int, c_int, mode_c)
y1 = get_pixel(<double*>img.data, rows, columns,
r_int + 1, c_int, mode_c)
y2 = get_pixel(<double*>img.data, rows, columns,
r_int + 1, c_int + 1, mode_c)
y3 = get_pixel(<double*>img.data, rows, columns,
r_int, c_int + 1, mode_c)
out[tfr, tfc] = \
(1 - t) * (1 - u) * y0 + \
t * (1 - u) * y1 + \
t * u * y2 + (1 - t) * u * y3;
out[tfr, tfc] = bilinear_interpolation(<double*>img.data, rows, cols, r, c, mode_c)
return out
+11 -12
View File
@@ -2,7 +2,7 @@ from numpy.testing import assert_array_almost_equal, run_module_suite
import numpy as np
from skimage.transform import (warp, homography, fast_homography,
SimilarityTransform)
SimilarityTransform, ProjectiveTransform)
from skimage import transform as tf, data, img_as_float
from skimage.color import rgb2gray
@@ -34,7 +34,7 @@ def test_homography():
def test_fast_homography():
img = rgb2gray(data.lena()).astype(np.uint8)
img = rgb2gray(data.lena())
img = img[:, :100]
theta = np.deg2rad(30)
@@ -49,20 +49,19 @@ def test_fast_homography():
H[:2, 2] = [tx, ty]
for mode in ('constant', 'mirror', 'wrap'):
p0 = homography(img, H, mode=mode, order=1)
p0 = warp(img, ProjectiveTransform(H).inverse, mode=mode, order=1)
p1 = fast_homography(img, H, mode=mode)
p1 = np.round(p1)
## import matplotlib.pyplot as plt
## f, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4)
## ax0.imshow(img)
## ax1.imshow(p0, cmap=plt.cm.gray)
## ax2.imshow(p1, cmap=plt.cm.gray)
## ax3.imshow(np.abs(p0 - p1), cmap=plt.cm.gray)
## plt.show()
# import matplotlib.pyplot as plt
# f, (ax0, ax1, ax2, ax3) = plt.subplots(1, 4)
# ax0.imshow(img)
# ax1.imshow(p0, cmap=plt.cm.gray)
# ax2.imshow(p1, cmap=plt.cm.gray)
# ax3.imshow(np.abs(p0 - p1), cmap=plt.cm.gray)
# plt.show()
d = np.mean(np.abs(p0 - p1))
assert d < 0.2
assert d < 0.001
def test_swirl():