Created
October 16, 2012 17:05
-
-
Save pmarshwx/3900591 to your computer and use it in GitHub Desktop.
Rotated, Anisotropic Gaussian Kernel Density Estimation
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
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 |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment