wassname/scikit-image
1
0
Fork 0
mirror of https://github.com/wassname/scikit-image.git synced 2026-08-05 13:21:12 +08:00
Code Issues Packages Projects Releases Wiki Activity

Cythonizing the for loops

Browse Source
This commit is contained in:
Ankit Agrawal
2013-08-15 14:59:54 +05:30
parent d15741c58c
commit 7c32dc0c15
2 changed files with 70 additions and 42 deletions
Download Patch File Download Diff File Expand all files Collapse all files
skimage/feature/censure.py
+10 -42
View File
@@ -5,11 +5,11 @@ from ..transform import integral_image
from ..feature.corner import _compute_auto_correlation
from ..util import img_as_float
from .censure_cy import _censure_dob_loop
from .censure_cy import _censure_dob_loop, _slanted_integral_image, _censure_octagon_loop
from time import time
def _get_filtered_image(image, no_of_scales, mode):
# TODO : Implement the STAR and Octagon mode
# TODO : Implement the STAR mode
if mode == 'DoB':
scales = np.zeros((image.shape[0], image.shape[1], no_of_scales))
for i in range(no_of_scales):
@@ -22,6 +22,7 @@ def _get_filtered_image(image, no_of_scales, mode):
scales[:, :, i] = filtered_image
return scales
elif mode == 'Octagon':
# TODO : Decide the shapes of Octagon filters for scales > 7
outer_shape = [(5, 2), (5, 3), (7, 3), (9, 4), (9, 7), (13, 7), (15, 10)]
inner_shape = [(3, 0), (3, 1), (3, 2), (5, 2), (5, 3), (5, 4), (5, 5)]
scales = np.zeros((image.shape[0], image.shape[1], no_of_scales))
@@ -41,48 +42,13 @@ def _get_filtered_image(image, no_of_scales, mode):
inner_pixels = (mi + 2 * ni)**2 - 2 * ni * (ni + 1)
outer_wt = 1.0 / (outer_pixels - inner_pixels)
inner_wt = 1.0 / inner_pixels
o_m = (mo - 1) / 2
i_m = (mi - 1) / 2
o_set = o_m + no
i_set = i_m + ni
# Outsource to Cython
for i in range(o_set + 1, image.shape[0] - o_set - 1):
for j in range(o_set + 1, image.shape[1] - o_set - 1):
outer = integral_img1[i + o_set, j + o_m] - integral_img1[i + o_m, j + o_set] - integral_img[i + o_set, j - o_m] + integral_img[i + o_m, j - o_m]
outer += integral_img[i + (mo - 3) / 2, j + (mo - 3) / 2] - integral_img[i - o_m, j + (mo - 3) / 2] - integral_img[i + (mo - 3) / 2, j - o_m] + integral_img[i - o_m, j - o_m]
outer += integral_img4[i + o_m, j - o_set] - integral_img4[i + o_set, j - o_m] - integral_img[i - o_m, j - (mo + 1) / 2] + integral_img[i - o_m, j - o_set - 1]
outer += integral_img2[i - o_set, j - o_m] - integral_img2[i - o_m, j - o_set] - integral_img[i - (mo + 1) / 2, -1] - integral_img[i - o_set - 1, j + (mo - 3) / 2] + integral_img[i - (mo + 1) / 2, j + (mo - 3) / 2] + integral_img[i - o_set - 1, -1]
outer += integral_img3[i - o_m, j + o_set] - integral_img3[i - o_set, j + o_m] - integral_img[-1, j + o_set + 1] - integral_img[i + (mo - 3) / 2, j + o_m] + integral_img[-1, j + o_m] + integral_img[i + (mo - 3) / 2, j + o_set + 1]
inner = integral_img1[i + i_set, j + i_m] - integral_img1[i + i_m, j + i_set] - integral_img[i + i_set, j - i_m] + integral_img[i + i_m, j - i_m]
inner += integral_img[i + (mi - 3) / 2, j + (mi - 3) / 2] - integral_img[i - i_m, j + (mi - 3) / 2] - integral_img[i + (mi - 3) / 2, j - i_m] + integral_img[i - i_m, j - i_m]
inner += integral_img4[i + i_m, j - i_set] - integral_img4[i + i_set, j - i_m] - integral_img[i - i_m, j - (mi + 1) / 2] + integral_img[i - i_m, j - i_set - 1]
inner += integral_img2[i - i_set, j - i_m] - integral_img2[i - i_m, j - i_set] - integral_img[i - (mi + 1) / 2, -1] - integral_img[i - i_set - 1, j + (mi - 3) / 2] + integral_img[i - (mi + 1) / 2, j + (mi - 3) / 2] + integral_img[i - i_set - 1, -1]
inner += integral_img3[i - i_m, j + i_set] - integral_img3[i - i_set, j + i_m] - integral_img[-1, j + i_set + 1] - integral_img[i + (mi - 3) / 2, j + i_m] + integral_img[-1, j + i_m] + integral_img[i + (mi - 3) / 2, j + i_set + 1]
_censure_octagon_loop(image, integral_img, integral_img1, integral_img2, integral_img3, integral_img4, filtered_image, outer_wt, inner_wt, mo, no, mi, ni)
filtered_image[i, j] = outer_wt * outer - (outer_wt + inner_wt) * inner
scales[:, :, k] = filtered_image
return scales
# Outsource to Cython
def _slanted_integral_image(image):
flipped_lr = np.fliplr(image)
left_sum = np.zeros(image.shape[0])
for i in range(image.shape[1] - image.shape[0], image.shape[1]):
left_sum[image.shape[1] - 1 - i] = np.sum(flipped_lr.diagonal(i))
left_sum = left_sum.cumsum(0)
right_sum = np.sum(image, 1).cumsum(0)
image[:, 0] = left_sum
image[:, -1] = right_sum
integral_img = np.zeros((image.shape[0] + 1, image.shape[1]))
integral_img[1:, :] = image
for i in range(1, integral_img.shape[0]):
for j in range(1, integral_img.shape[1] - 1):
integral_img[i, j] += integral_img[i, j - 1] + integral_img[i - 1, j + 1] - integral_img[i - 1, j]
return integral_img[1:, :integral_img.shape[1]]
def _slanted_integral_image_modes(img, mode=1):
if mode == 1:
image = np.copy(img)
@@ -135,22 +101,24 @@ def censure_keypoints(image, no_of_scales=7, mode='DoB', threshold=0.03, rpc_thr
image = np.ascontiguousarray(image)
# Generating all the scales
start = time()
scales = np.zeros((image.shape[0], image.shape[1], no_of_scales))
scales = _get_filtered_image(image, no_of_scales, mode)
print time() - start
# Suppressing points that are neither minima or maxima in their 3 x 3 x 3
# neighbourhood to zero
minimas = (minimum_filter(scales, (3, 3, 3)) == scales).astype(int) * scales
maximas = (maximum_filter(scales, (3, 3, 3)) == scales).astype(int) * scales
print time() - start
# Suppressing minimas and maximas weaker than threshold
minimas[np.abs(minimas) < threshold] = 0
maximas[np.abs(maximas) < threshold] = 0
response = maximas + minimas
print time() - start
for i in range(1, no_of_scales - 1):
response[:, :, i] = _suppress_line(response[:, :, i], (1 + i / 3.0), rpc_threshold)
print time() - start
# Returning keypoints with its scale
keypoints = np.transpose(np.nonzero(response[:, :, 1:no_of_scales])) + [0, 0, 1]
return keypoints
skimage/feature/censure_cy.pyx
+60
View File
@@ -4,6 +4,7 @@
#cython: wraparound=False
cimport numpy as cnp
import numpy as np
def _censure_dob_loop(double[:, ::1] image, cnp.int16_t n,
@@ -19,3 +20,62 @@ def _censure_dob_loop(double[:, ::1] image, cnp.int16_t n,
inner = integral_img[i + n, j + n] + integral_img[i - n - 1, j - n - 1] - integral_img[i + n, j - n - 1] - integral_img[i - n - 1, j + n]
outer = integral_img[i + 2 * n, j + 2 * n] + integral_img[i - 2 * n - 1, j - 2 * n - 1] - integral_img[i + 2 * n, j - 2 * n - 1] - integral_img[i - 2 * n - 1, j + 2 * n]
filtered_image[i, j] = outer_wt * outer - (inner_wt + outer_wt) * inner
def _slanted_integral_image(double[:, :] image):
cdef Py_ssize_t i, j
cdef double[:, :] flipped_lr = np.fliplr(image)
cdef double[:] left_sum = np.zeros(image.shape[0], dtype=np.float)
cdef double[:] right_sum
cdef double[:, :] integral_img = np.zeros((image.shape[0] + 1, image.shape[1]), dtype=np.float)
flipped_lr = np.fliplr(image)
for i in range(image.shape[1] - image.shape[0], image.shape[1]):
left_sum[image.shape[1] - 1 - i] = np.sum(flipped_lr.diagonal(i))
left_sum = left_sum.cumsum(0)
right_sum = np.sum(image, 1).cumsum(0)
image[:, 0] = left_sum
image[:, -1] = right_sum
integral_img[1:, :] = image
for i in range(1, integral_img.shape[0]):
for j in range(1, integral_img.shape[1] - 1):
integral_img[i, j] += integral_img[i, j - 1] + integral_img[i - 1, j + 1] - integral_img[i - 1, j]
return integral_img[1:, :integral_img.shape[1]]
def _censure_octagon_loop(double[:, ::1] image, double[:, ::1] integral_img,
double[:, ::1] integral_img1,
double[:, ::1] integral_img2,
double[:, ::1] integral_img3,
double[:, ::1] integral_img4,
double[:, ::1] filtered_image,
double outer_wt, double inner_wt,
int mo, int no, int mi, int ni):
cdef Py_ssize_t i, j, o_m, i_m, o_set, i_set
o_m = (mo - 1) / 2
i_m = (mi - 1) / 2
o_set = o_m + no
i_set = i_m + ni
for i in range(o_set + 1, image.shape[0] - o_set - 1):
for j in range(o_set + 1, image.shape[1] - o_set - 1):
outer = integral_img1[i + o_set, j + o_m] - integral_img1[i + o_m, j + o_set] - integral_img[i + o_set, j - o_m] + integral_img[i + o_m, j - o_m]
outer += integral_img[i + o_m - 1, j + o_m - 1] - integral_img[i - o_m, j + o_m - 1] - integral_img[i + o_m - 1, j - o_m] + integral_img[i - o_m, j - o_m]
outer += integral_img4[i + o_m, j - o_set] - integral_img4[i + o_set, j - o_m] - integral_img[i - o_m, j - o_m + 1] + integral_img[i - o_m, j - o_set - 1]
outer += integral_img2[i - o_set, j - o_m] - integral_img2[i - o_m, j - o_set] - integral_img[i - o_m + 1, -1] - integral_img[i - o_set - 1, j + o_m - 1] + integral_img[i - o_m + 1, j + o_m - 1] + integral_img[i - o_set - 1, -1]
outer += integral_img3[i - o_m, j + o_set] - integral_img3[i - o_set, j + o_m] - integral_img[-1, j + o_set + 1] - integral_img[i + o_m - 1, j + o_m] + integral_img[-1, j + o_m] + integral_img[i + o_m - 1, j + o_set + 1]
inner = integral_img1[i + i_set, j + i_m] - integral_img1[i + i_m, j + i_set] - integral_img[i + i_set, j - i_m] + integral_img[i + i_m, j - i_m]
inner += integral_img[i + i_m - 1, j + i_m - 1] - integral_img[i - i_m, j + i_m - 1] - integral_img[i + i_m - 1, j - i_m] + integral_img[i - i_m, j - i_m]
inner += integral_img4[i + i_m, j - i_set] - integral_img4[i + i_set, j - i_m] - integral_img[i - i_m, j - i_m + 1] + integral_img[i - i_m, j - i_set - 1]
inner += integral_img2[i - i_set, j - i_m] - integral_img2[i - i_m, j - i_set] - integral_img[i - i_m + 1, -1] - integral_img[i - i_set - 1, j + i_m - 1] + integral_img[i - i_m + 1, j + i_m - 1] + integral_img[i - i_set - 1, -1]
inner += integral_img3[i - i_m, j + i_set] - integral_img3[i - i_set, j + i_m] - integral_img[-1, j + i_set + 1] - integral_img[i + i_m - 1, j + i_m] + integral_img[-1, j + i_m] + integral_img[i + i_m - 1, j + i_set + 1]
filtered_image[i, j] = outer_wt * outer - (outer_wt + inner_wt) * inner