public

Rotated, Anisotropic Gaussian Kernel Density Estimation

  • Download Gist
rotated_anisotropic_gaussian_kde.pyx
Cython
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101
cimport cython
import numpy as np
cimport numpy as np
 
cdef extern from 'math.h':
float exp(float x)
float cos(float x)
float sin(float x)
float fabs(float x)
 
DTYPE64 = np.float64
ctypedef np.float64_t DTYPE64_t
 
 
@cython.boundscheck(False)
@cython.cdivision(True)
cpdef anisotropic_gauss(np.ndarray[DTYPE64_t, ndim=2] data,
float sigx, float sigy, float rot,
int h, int k, float dx, float factor,
sig_as_grid_points=False):
 
cdef unsigned int vlength = data.shape[0]
cdef unsigned int ulength = data.shape[1]
cdef unsigned int a, b, nxy, pwdimlength
cdef int nxyn, iw, ie, js, jn, nw
cdef float sigx_sq, sigy_sq, rad, amp
cdef float sintheta, sintheta_sq, costheta, costheta_sq
cdef float A, B, C, D, E, F, ellipse, X, Y
cdef Py_ssize_t i, j, ii, jj, x, y
 
cdef np.ndarray[DTYPE64_t, ndim=2] frc_data = np.zeros([vlength, ulength], dtype=DTYPE64)
cdef float PI=3.141592653589793
 
if not sig_as_grid_points:
sigx = float(sigx / dx)
sigy = float(sigy / dx)
 
sigx_sq = sigx*sigx
sigy_sq = sigy*sigy
 
# Set up X-axis
a = int(factor * sigx)
a_sq = float(a * a)
 
# Set up Y-axis
b = int(factor * sigy)
b_sq = float(b * b)
 
# Set weight-loop to size of maximum length
if a >= b:
nxy = a
else:
nxy = b
 
nxyn = -1 * nxy
pwdimlength = 2*nxy+1
 
# Rotation Matrix Values
rad = rot * PI / 180. * -1. # Not entirely sure why the -1, but it makes it work
sintheta = sin(rad)
sintheta_sq = sintheta**2
costheta = cos(rad)
costheta_sq = costheta**2
 
# Constants used in transformation
A = (costheta_sq / a_sq) + (sintheta_sq / b_sq)
B = -2 * costheta * sintheta * ((1 / a_sq) - (1 / b_sq))
C = (sintheta_sq / a_sq) + (costheta_sq / b_sq)
# D = (2 * A * h ) + (B * k)
# E = (2 * C * k) + (B * h)
# F = (A * h * h) + (B * h * k) + (C * k * k)
 
cdef np.ndarray[DTYPE64_t, ndim=1] partweight = np.zeros([pwdimlength**2], dtype=DTYPE64)
 
nw = -1
for y in range(nxyn, nxy+1):
for x in range(nxyn, nxy+1):
nw += 1
ellipse = (A*x*x) + (B*x*y) + (C*y*y) # - D*X) - (E*Y) + F
if ellipse <= 1:
X = (x)*costheta - (y)*sintheta
Y = (x)*sintheta + (y)*costheta
partweight[nw] = exp(-0.5*((X*X / sigx_sq) + (Y*Y / sigy_sq)))
 
for j in range(0, vlength):
for i in range(0, ulength):
if data[j,i] != 0:
amp = data[j,i] / (2*PI*sigx*sigy)
iw=i-nxy
ie=i+nxy
js=j-nxy
jn=j+nxy
nw = -1
for jj in range(js+k, jn+k+1):
for ii in range(iw+h, ie+h+1):
nw += 1
if jj < 0 or jj >= vlength or ii < 0 or ii >= ulength:
continue
frc_data[jj,ii] = frc_data[jj,ii] + amp*partweight[nw]
 
return frc_data

Please sign in to comment on this gist.

Something went wrong with that request. Please try again.