Instantly share code, notes, and snippets.

# pmarshwx/rotated_anisotropic_gaussian_kde.pyx Created Oct 16, 2012

Rotated, Anisotropic Gaussian Kernel Density Estimation
 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
to join this conversation on GitHub. Already have an account? Sign in to comment