New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

WIP: multi-threaded on affine registration #1571

Closed
wants to merge 4 commits into
base: master
from
Closed
Diff settings

Always

Just for now

Copy path View file
@@ -8,6 +8,8 @@ cimport numpy as cnp
cimport cython
from .fused_types cimport floating, number

from cython.parallel import prange, parallel
from libc.stdlib cimport abort, malloc, free

cdef extern from "dpy_math.h" nogil:
double floor(double)
@@ -3104,25 +3106,32 @@ def _gradient_3d(floating[:, :, :] img, double[:, :] img_world2grid,
int nslices = out.shape[0]
int nrows = out.shape[1]
int ncols = out.shape[2]
int i, j, k, in_flag
double tmp
double[:] x = np.empty(shape=(3,), dtype=np.float64)
double[:] dx = np.empty(shape=(3,), dtype=np.float64)
int i, j, k, in_flag, p
double tmp, *x, *dx, *q
double[:] h = np.empty(shape=(3,), dtype=np.float64)
double[:] q = np.empty(shape=(3,), dtype=np.float64)
with nogil:
with nogil, parallel():
h[0] = 0.5 * img_spacing[0]
h[1] = 0.5 * img_spacing[1]
h[2] = 0.5 * img_spacing[2]
for k in range(nslices):
x = <double *>malloc(sizeof(double) * 3)
if x == NULL:
abort()
dx = <double *>malloc(sizeof(double) * 3)
if dx == NULL:
abort()
q = <double *>malloc(sizeof(double) * 3)
if q == NULL:
abort()
for k in prange(nslices, schedule='guided'):
for i in range(nrows):
for j in range(ncols):
inside[k, i, j] = 1
# Compute coordinates of index (k, i, j) in physical space
x[0] = _apply_affine_3d_x0(k, i, j, 1, out_grid2world)
x[1] = _apply_affine_3d_x1(k, i, j, 1, out_grid2world)
x[2] = _apply_affine_3d_x2(k, i, j, 1, out_grid2world)
dx[:] = x[:]
for p in range(3):
dx[p] = x[p]
for p in range(3):
# Compute coordinates of point dx on img's grid
dx[p] = x[p] - h[p]
@@ -3157,6 +3166,9 @@ def _gradient_3d(floating[:, :, :] img, double[:, :] img_world2grid,
continue
out[k, i, j, p] = (out[k, i, j, p] - tmp) / img_spacing[p]
dx[p] = x[p]
free(x)
free(dx)
free(q)


def _sparse_gradient_3d(floating[:, :, :] img,
@@ -3194,16 +3206,20 @@ def _sparse_gradient_3d(floating[:, :, :] img,
"""
cdef:
int n = sample_points.shape[0]
int i, in_flag
double tmp
double[:] dx = np.empty(shape=(3,), dtype=np.float64)
int i, in_flag, p
double tmp, *dx, *q
double[:] h = np.empty(shape=(3,), dtype=np.float64)
double[:] q = np.empty(shape=(3,), dtype=np.float64)
with nogil:
with nogil, parallel():
h[0] = 0.5 * img_spacing[0]
h[1] = 0.5 * img_spacing[1]
h[2] = 0.5 * img_spacing[2]
for i in range(n):
dx = <double *>malloc(sizeof(double) * 3)
if dx == NULL:
abort()
q = <double *>malloc(sizeof(double) * 3)
if q == NULL:
abort()
for i in prange(n, schedule='guided'):
inside[i] = 1
dx[0] = sample_points[i, 0]
dx[1] = sample_points[i, 1]
@@ -3242,6 +3258,8 @@ def _sparse_gradient_3d(floating[:, :, :] img,
continue
out[i, p] = (out[i, p] - tmp) / img_spacing[p]
dx[p] = sample_points[i, p]
free(dx)
free(q)


def _gradient_2d(floating[:, :] img, double[:, :] img_world2grid,
@@ -3283,7 +3301,7 @@ def _gradient_2d(floating[:, :] img, double[:, :] img_world2grid,
cdef:
int nrows = out.shape[0]
int ncols = out.shape[1]
int i, j, k, in_flag
int i, j, k, in_flag, p
double tmp
double[:] x = np.empty(shape=(2,), dtype=np.float64)
double[:] dx = np.empty(shape=(2,), dtype=np.float64)
@@ -3363,7 +3381,7 @@ def _sparse_gradient_2d(floating[:, :] img, double[:, :] img_world2grid,
"""
cdef:
int n = sample_points.shape[0]
int i, in_flag
int i, in_flag, p
double tmp
double[:] dx = np.empty(shape=(2,), dtype=np.float64)
double[:] h = np.empty(shape=(2,), dtype=np.float64)
ProTip! Use n and p to navigate between commits in a pull request.