From d7979cc2759d139d2e3c7ab084fa6887b3e9cb4d Mon Sep 17 00:00:00 2001
From: "Joern Ungermann (IEK-7 FZ-Juelich)" <j.ungermann@fz-juelich.de>
Date: Thu, 28 Nov 2013 15:12:33 +0100
Subject: [PATCH] Added numerical core routines for Jacobian multiply.
---
pyramid/numcore/phase_mag_real.pyx | 63 +++++++++++++++++++++++++++++-
pyramid/phasemapper.py | 16 ++++++++
2 files changed, 77 insertions(+), 2 deletions(-)
diff --git a/pyramid/numcore/phase_mag_real.pyx b/pyramid/numcore/phase_mag_real.pyx
index d679933..a0403f5 100644
--- a/pyramid/numcore/phase_mag_real.pyx
+++ b/pyramid/numcore/phase_mag_real.pyx
@@ -92,8 +92,8 @@ def get_jacobi_core(
@cython.wraparound(False)
def get_jacobi_core(
unsigned int v_dim, unsigned int u_dim,
- np.ndarray v_phi, np.ndarray u_phi,
- np.ndarray jacobi):
+ double[:, :] v_phi, double[:, :] u_phi,
+ double[:, :] jacobi):
'''DOCSTRING!'''
# TODO: Docstring!!!
cdef unsigned int i, j, p, q, p_min, p_max, q_min, q_max, u_column, v_column, row
@@ -105,3 +105,62 @@ def get_jacobi_core(
# v_dim*u_dim columns for the v-component (note the minus!):
jacobi[:, u_dim*v_dim+i+u_dim*j] = \
-np.reshape(v_phi[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i], -1)
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def multiply_jacobi_core(
+ unsigned int v_dim, unsigned int u_dim,
+ double[:, :] v_phi, double[:, :] u_phi,
+ double[:] vector,
+ double[:] result):
+
+ cdef unsigned int s, i, j, u_min, u_max, v_min, v_max, ri, u, v, siz
+ cdef double v0, v1
+ siz = u_dim * v_dim
+
+ s = 0
+ for i in range(u_dim):
+ for j in range(v_dim):
+ u_min = (u_dim - 1) - i
+ v_min = (v_dim - 1) - j
+
+ ri = 0
+ for v in range(v_min, v_min + v_dim):
+ for u in range(u_min, u_min + u_dim):
+ result[ri] += vector[s] * u_phi[v, u]
+ result[ri] -= vector[s + siz] * v_phi[v, u]
+ ri += 1
+ s += 1
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def multiply_jacobi_core2(
+ int v_dim, int u_dim,
+ double[:, :] v_phi, double[:, :] u_phi,
+ double[:] vector,
+ double[:] result):
+
+ cdef int s1, s2, i, j, u_min, u_max, v_min, v_max, ri, u, v, siz, j_min, j_max, i_min, i_max
+
+ siz = u_dim * v_dim
+
+ for v in range(2 * v_dim - 1):
+ for u in range(2 * u_dim - 1):
+ i_min = max(0, (u_dim - 1) - u)
+ i_max = min(u_dim, ((2 * u_dim) - 1) - u)
+ for i in range(i_min, i_max):
+ s1 = i * u_dim
+ s2 = s1 + siz
+
+ j_min = max(0, (v_dim - 1) - v)
+ j_max = min(v_dim, ((2 * v_dim) - 1) - v)
+
+ u_min = (u_dim - 1) - i
+ v_min = (v_dim - 1) - j_min
+ ri = (v - ((v_dim - 1) - j_min)) * u_dim + (u - u_min)
+ for j in range(j_min, j_max):
+ result[ri] += vector[s1 + j] * u_phi[v, u]
+ result[ri] -= vector[s2 + j] * v_phi[v, u]
+ ri += u_dim
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index 8ddf90b..fdc41c0 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -99,6 +99,22 @@ class Kernel:
result += vector[s+size]*-self.v[v_min:v_max, u_min:u_max].reshape(-1) # v
return result
+ def multiply_jacobi_core(self, vector):
+ v_dim, u_dim = self.dim
+ size = v_dim * u_dim
+ result = np.zeros(size)
+ nc.multiply_jacobi_core(
+ v_dim, u_dim, self.v, self.u, vector, result)
+ return result
+
+ def multiply_jacobi_core2(self, vector):
+ v_dim, u_dim = self.dim
+ size = v_dim * u_dim
+ result = np.zeros(size)
+ nc.multiply_jacobi_core2(
+ v_dim, u_dim, self.v, self.u, vector, result)
+ return result
+
def multiply_jacobi_T(self, vector):
# TODO: Docstring!!!
# vector: v_dim*u_dim elements for u_mag and v_dim*u_dim elements for v_mag
--
GitLab