hog cython
pip安装cython之后，将下面代码写入hogtest2.pyx文件（我通过改文件后缀新建）
import numpy as np
from PIL import Image
cimport numpy as np
cimport cython
cdef extern from "math.h":
double sqrt(double i)
double fabs(double i)
double floor(double i)
cpdef hog(im, int sbin = 8):
cdef np.ndarray[np.double_t, ndim=3] data, feat
cdef np.ndarray[np.double_t, ndim=1] hist, norm
cdef int blocks0, blocks1
cdef int out0, out1, out2
cdef int visible0, visible1
cdef double dy, dx, v
cdef double dy2, dx2, v2
cdef double dy3, dx3, v3
cdef double best_dot, dot
cdef int best_o
cdef double xp, yp, vx0, vy0, vx1, vy1
cdef int ixp, iyp
cdef double n1, n2, n3, n4, t1, t2, t3, t4, h1, h2, h3, h4
cdef int p
cdef np.ndarray[np.double_t, ndim=1] uu
uu = np.array([ 1.0000, 0.9397, 0.7660, 0.500, 0.1736,
-0.1736, -0.5000, -0.7660, -0.9397])
cdef np.ndarray[np.double_t, ndim=1] vv
vv = np.array([0.0000, 0.3420, 0.6428, 0.8660, 0.9848,
0.9848, 0.8660, 0.6428, 0.3420])
cdef double eps = 0.0001 # to avoid division by 0
cdef unsigned int cc0, cc1, cc2
cdef int x, y, o, q
cdef int dstptr, srcptr
width, height = im.size
blocks0 = height / sbin
blocks1 = width / sbin
out0 = blocks0 - 2
out1 = blocks1 - 2
out2 = 9 + 4
visible0 = blocks0 * sbin
visible1 = blocks1 * sbin
data = np.asarray(im, dtype=np.double)
cc0 = <unsigned int>(0)
cc1 = <unsigned int>(1)
cc2 = <unsigned int>(2)
hist = np.zeros(shape=(blocks0 * blocks1 * 9), dtype=np.double)
norm = np.zeros(shape=(blocks0 * blocks1), dtype=np.double)
feat = np.zeros(shape=(out0, out1, out2), dtype=np.double)
for x from 1 <= x < visible1 - 1:
for y from 1 <= y < visible0 - 1:
dy = data[y + 1, x, cc0] - data[y - 1, x, cc0]
dx = data[y, x + 1, cc0] - data[y, x - 1, cc0]
v = dx * dx + dy * dy
dy2 = data[y + 1, x, cc1] - data[y - 1, x, cc1]
dx2 = data[y, x + 1, cc1] - data[y, x - 1, cc1]
v2 = dx2 * dx2 + dy2 * dy2
dy3 = data[y + 1, x, cc2] - data[y - 1, x, cc2]
dx3 = data[y, x + 1, cc2] - data[y, x - 1, cc2]
v3 = dx3 * dx3 + dy3 * dy3
if v2 > v: # pick channel with strongest gradient
v = v2
dx = dx2
dy = dy2
if v3 > v:
v = v3
dx = dx3
dy = dy3
# snap to one of 9 orientations
best_dot = 0.
best_o = 0
for o from 0 <= o < 9:
dot = fabs(uu[o] * dx + vv[o] * dy)
if dot > best_dot:
best_dot = dot
best_o = o
# add to 4 histograms around pixel using linear interpolation
xp = (<double>(x) + 0.5) / <double>(sbin) - 0.5
yp = (<double>(y) + 0.5) / <double>(sbin) - 0.5
ixp = <int>floor(xp)
iyp = <int>floor(yp)
vx0 = xp - ixp
vy0 = yp - iyp
vx1 = 1.0 - vx0
vy1 = 1.0 - vy0
v = sqrt(v)
if ixp >= 0 and iyp >= 0:
hist[ixp * blocks0 + iyp + best_o*blocks0*blocks1] += vx1 * vy1 * v
if ixp + 1 < blocks1 and iyp >= 0:
hist[(ixp + 1) * blocks0 + iyp + best_o*blocks0*blocks1] += vx0 * vy1 * v
if ixp >= 0 and iyp + 1 < blocks0:
hist[ixp * blocks0 + (iyp + 1) + best_o*blocks0*blocks1] += vx1 * vy0 * v
if ixp + 1 < blocks1 and iyp + 1 < blocks0:
hist[(ixp + 1) * blocks0 + (iyp + 1) + best_o * blocks0 * blocks1] += vx0 * vy0 * v
# compute energy in each block by summing over orientations
for o from 0 <= o < 9:
for q from 0 <= q < blocks0 * blocks1:
norm[q] += hist[o * blocks0 * blocks1 + q] * hist[o * blocks0 * blocks1 + q]
# compute normalized values
for x from 0 <= x < out1:
for y from 0 <= y < out0:
p = (x+1) * blocks0 + y + 1
n1 = 1.0 / sqrt(norm[p] + norm[p+1] + norm[p+blocks0] + norm[p+blocks0+1] + eps)
p = (x+1) * blocks0 + y
n2 = 1.0 / sqrt(norm[p] + norm[p+1] + norm[p+blocks0] + norm[p+blocks0+1] + eps)
p = x * blocks0 + y + 1
n3 = 1.0 / sqrt(norm[p] + norm[p+1] + norm[p+blocks0] + norm[p+blocks0+1] + eps)
p = x * blocks0 + y
n4 = 1.0 / sqrt(norm[p] + norm[p+1] + norm[p+blocks0] + norm[p+blocks0+1] + eps)
t1 = 0
t2 = 0
t3 = 0
t4 = 0
srcptr = (x+1) * blocks0 + y + 1
for o from 0 <= o < 9:
h1 = hist[srcptr] * n1
h2 = hist[srcptr] * n2
h3 = hist[srcptr] * n3
h4 = hist[srcptr] * n4
# for some reason, gcc will not automatically inline
# the min function here, so we just do it ourselves
# for impressive speedups
if h1 > 0.2:
h1 = 0.2
if h2 > 0.2:
h2 = 0.2
if h3 > 0.2:
h3 = 0.2
if h4 > 0.2:
h4 = 0.2
feat[y, x, o] = 0.5 * (h1 + h2 + h3 + h4)
t1 += h1
t2 += h2
t3 += h3
t4 += h4
srcptr += blocks0 * blocks1
feat[y, x, 9] = 0.2357 * t1
dstptr += out0 * out1
feat[y, x, 10] = 0.2357 * t2
dstptr += out0 * out1
feat[y, x, 11] = 0.2357 * t3
dstptr += out0 * out1
feat[y, x, 12] = 0.2357 * t4
return feat
cpdef hogpad(np.ndarray[np.double_t, ndim=3] hog):
cdef np.ndarray[np.double_t, ndim=3] out
cdef int i, j, k
cdef int w = hog.shape[0], h = hog.shape[1], z = hog.shape[2]
out = np.zeros((w + 2, h + 2, z))
for i in range(w):
for j in range(h):
for k in range(z):
out[i+1, j+1, k] = hog[i, j, k]
return out
cpdef rgbhist(im, int binsize = 8):
"""
Computes an RGB color histogram with a binsize.
"""
cdef int w = im.size[0], h = im.size[1]
cdef np.ndarray[np.uint8_t, ndim=3] data = np.asarray(im)
cdef np.ndarray[np.double_t, ndim=1] hist
hist = np.zeros(binsize * binsize * binsize)
for i from 0 <= i < w:
for j from 0 <= j < h:
bin = (<int>data[j,i,0]) / (256/binsize)
bin += (<int>data[j,i,1]) / (256/binsize) * binsize
bin += (<int>data[j,i,2]) / (256/binsize) * binsize * binsize
hist[bin] += 1
return hist
cpdef rgbmean(im):
"""
Computes mean and covariances of RGB colors.
"""
cdef int w = im.size[0], h = im.size[1]
cdef double r, g, b
cdef np.ndarray[np.uint8_t, ndim=3] data = np.asarray(im)
cdef np.ndarray[np.double_t, ndim=1] out = np.zeros(9)
for i in range(w):
for j in range(h):
r = data[j, i, 0] / 255.
g = data[j, i, 1] / 255.
b = data[j, i, 2] / 255.
out[0] += r
out[1] += g
out[2] += b
out[3] += r * r
out[4] += r * g
out[5] += r * b
out[6] += g * g
out[7] += g * b
out[8] += b * b
return out / (w * h)
然后，新建setup.py
from distutils.core import setup
from distutils.extension import Extension
from Cython.Build import cythonize
import numpy
setup(name='hog2',
ext_modules=cythonize("hogtest2.pyx"),
include_dirs=[numpy.get_include()]
)
cmd 运行 python setup.py build
继续运行python setup.py install
将运行后出现的build文件夹中的pyd后缀文件放到python的LIB\sit-packages里，或者放在待会要跑的py文件同目录下
我新建了一个test.py
from PIL import Image
import hogtest2
img = Image.open("cat.jpg")
hogtest2.hog(img)
a=hogtest2.rgbmean(img)
print(a)
是可用的
