Skip to content

Commit

Permalink
Merge branch 'dev/0.2.3a0'
Browse files Browse the repository at this point in the history
  • Loading branch information
Malcolm White committed May 22, 2020
2 parents a0b2a15 + 304f33d commit e407e0c
Show file tree
Hide file tree
Showing 5 changed files with 210 additions and 48 deletions.
9 changes: 9 additions & 0 deletions pykonal/__init__.py
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,12 @@
from .__version__ import (
__major_version__,
__minor_version__,
__patch__,
__release__,
__version_tuple__,
__version__
)

from .solver import EikonalSolver

from . import constants
Expand Down
12 changes: 11 additions & 1 deletion pykonal/__version__.py
View file
Original file line number Diff line number Diff line change
@@ -1 +1,11 @@
__version__ = '0.2.1b0'
__major_version__ = 0
__minor_version__ = 2
__patch__ = 3
__release__ = "a0"
__version_tuple__ = (
__major_version__,
__minor_version__,
__patch__
)
__version_number__ = ".".join([str(v) for v in __version_tuple__])
__version__ = f"{__version_number__}{__release__}"
36 changes: 26 additions & 10 deletions pykonal/fields.pxd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,13 +4,14 @@ from . cimport constants


cdef class Field3D(object):
cdef str cy_coord_sys
cdef constants.REAL_t[3] cy_max_coords
cdef constants.REAL_t[3] cy_min_coords
cdef constants.UINT_t[3] cy_npts
cdef constants.REAL_t[3] cy_node_intervals
cdef constants.BOOL_t[3] cy_iax_isnull
cdef constants.BOOL_t[3] cy_iax_isperiodic
cdef str cy_coord_sys
cdef constants.BOOL_t[3] cy_iax_isnull
cdef constants.BOOL_t[3] cy_iax_isperiodic
cdef constants.REAL_t[3] cy_max_coords
cdef constants.REAL_t[3] cy_min_coords
cdef constants.REAL_t[:,:,:,:] cy_norm
cdef constants.UINT_t[3] cy_npts
cdef constants.REAL_t[3] cy_node_intervals

cdef constants.BOOL_t _update_max_coords(Field3D self)
cdef constants.BOOL_t _update_iax_isnull(Field3D self)
Expand All @@ -20,16 +21,31 @@ cdef class Field3D(object):
cdef class ScalarField3D(Field3D):
cdef constants.REAL_t[:,:,:] cy_values

cpdef np.ndarray[constants.REAL_t, ndim=1] resample(ScalarField3D self, constants.REAL_t[:,:] points, constants.REAL_t null=*)
cpdef constants.REAL_t value(ScalarField3D self, constants.REAL_t[:] point, constants.REAL_t null=*)
cpdef np.ndarray[constants.REAL_t, ndim=1] resample(
ScalarField3D self,
constants.REAL_t[:,:] points,
constants.REAL_t null=*
)
cpdef np.ndarray[constants.REAL_t, ndim=2] trace_ray(
ScalarField3D self,
constants.REAL_t[:] end
)
cpdef constants.REAL_t value(
ScalarField3D self,
constants.REAL_t[:] point,
constants.REAL_t null=*
)
cpdef VectorField3D _gradient_of_cartesian(ScalarField3D self)
cpdef VectorField3D _gradient_of_spherical(ScalarField3D self)


cdef class VectorField3D(Field3D):
cdef constants.REAL_t[:,:,:,:] cy_values

cpdef np.ndarray[constants.REAL_t, ndim=1] value(VectorField3D self, constants.REAL_t[:] point)
cpdef np.ndarray[constants.REAL_t, ndim=1] value(
VectorField3D self,
constants.REAL_t[:] point
)


cpdef Field3D load(str path)
180 changes: 143 additions & 37 deletions pykonal/fields.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,10 @@ import numpy as np
from . import constants
from . import transformations

# Cython built-in imports.
from libc.math cimport sqrt, sin
from libcpp.vector cimport vector as cpp_vector
from libc.stdlib cimport malloc, free

# Third-party Cython imports
cimport numpy as np
Expand Down Expand Up @@ -76,41 +80,6 @@ cdef class Field3D(object):
)
return (arr)

@property
def node_intervals(self):
"""
[*Read/Write*, :class:`numpy.ndarray`\ (shape=(3,), dtype=numpy.float)]
Array of node intervals along each axis. This attribute must be
initialized by the user.
"""
return (np.asarray(self.cy_node_intervals))

@node_intervals.setter
def node_intervals(self, value):
value = np.asarray(value, dtype=constants.DTYPE_REAL)
if np.any(value) <= 0:
raise (ValueError("All node intervals must be > 0"))
self.cy_node_intervals = value
self._update_max_coords()
self._update_iax_isperiodic()


@property
def npts(self):
"""
[*Read/Write*, :class:`numpy.ndarray`\ (shape=(3,), dtype=numpy.int)]
Array specifying the number of nodes along each axis. This
attribute must be initialized by the user.
"""
return (np.asarray(self.cy_npts))

@npts.setter
def npts(self, value):
self.cy_npts = np.asarray(value, dtype=constants.DTYPE_UINT)
self._update_max_coords()
self._update_iax_isnull()
self._update_iax_isperiodic()

@property
def min_coords(self):
"""
Expand All @@ -128,7 +97,6 @@ cdef class Field3D(object):
self._update_max_coords()
self._update_iax_isperiodic()


@property
def max_coords(self):
"""
Expand All @@ -137,6 +105,23 @@ cdef class Field3D(object):
"""
return (np.asarray(self.cy_max_coords))

@property
def node_intervals(self):
"""
[*Read/Write*, :class:`numpy.ndarray`\ (shape=(3,), dtype=numpy.float)]
Array of node intervals along each axis. This attribute must be
initialized by the user.
"""
return (np.asarray(self.cy_node_intervals))

@node_intervals.setter
def node_intervals(self, value):
value = np.asarray(value, dtype=constants.DTYPE_REAL)
if np.any(value) <= 0:
raise (ValueError("All node intervals must be > 0"))
self.cy_node_intervals = value
self._update_max_coords()
self._update_iax_isperiodic()

@property
def nodes(self):
Expand All @@ -158,6 +143,52 @@ cdef class Field3D(object):
nodes = np.moveaxis(nodes, 0, -1)
return (nodes)

@property
def norm(self):
"""
[*Read-only*, :class:`numpy.ndarray`\ (shape=(N0,N1,N2,3), dtype=numpy.float)] 4D array of scaling
factors for gradient operator.
.. added:: 0.2.3a0
"""

try:
return (np.asarray(self.cy_norm))
except AttributeError:
norm = np.tile(
self.node_intervals,
np.append(self.npts, 1)
)
if self.coord_sys == "spherical":
norm[..., 1] *= self.nodes[..., 0]
norm[..., 2] *= self.nodes[..., 0]
norm[..., 2] *= np.sin(self.nodes[..., 1])
self.cy_norm = norm
return (np.asarray(self.cy_norm))

@property
def npts(self):
"""
[*Read/Write*, :class:`numpy.ndarray`\ (shape=(3,), dtype=numpy.int)]
Array specifying the number of nodes along each axis. This
attribute must be initialized by the user.
"""
return (np.asarray(self.cy_npts))

@npts.setter
def npts(self, value):
self.cy_npts = np.asarray(value, dtype=constants.DTYPE_UINT)
self._update_max_coords()
self._update_iax_isnull()
self._update_iax_isperiodic()

@property
def step_size(self):
"""
[*Read only*, :class:`float`] Step size used for ray tracing.
"""
return (self.norm[~np.isclose(self.norm, 0)].min() / 4)


cdef constants.BOOL_t _update_iax_isperiodic(Field3D self):
if self.cy_coord_sys == "spherical":
Expand Down Expand Up @@ -301,7 +332,82 @@ cdef class ScalarField3D(Field3D):
return (resampled)


cpdef constants.REAL_t value(ScalarField3D self, constants.REAL_t[:] point, constants.REAL_t null=np.nan):
cpdef np.ndarray[constants.REAL_t, ndim=2] trace_ray(
ScalarField3D self,
constants.REAL_t[:] end
):
"""
trace_ray(self, end)
Trace the ray ending at *end*.
This method traces the ray that ends at *end* in reverse
direction by taking small steps along the path of steepest
descent. The resulting ray is reversed before being returned,
so it is in the normal forward-time orientation.
:param end: Coordinates of the ray's end point.
:type end: numpy.ndarray(shape=(3,), dtype=numpy.float)
:return: The ray path ending at *end*.
:rtype: numpy.ndarray(shape=(N,3), dtype=numpy.float)
"""

cdef cpp_vector[constants.REAL_t *] ray
cdef constants.REAL_t norm, step_size, value, value_1back
cdef constants.REAL_t *point_new
cdef constants.REAL_t[3] gg, point, point_1back
cdef Py_ssize_t idx, jdx
cdef np.ndarray[constants.REAL_t, ndim=2] ray_np
cdef str coord_sys
cdef VectorField3D grad

coord_sys = self.coord_sys
step_size = self.step_size
grad = self.gradient

point_new = <constants.REAL_t *> malloc(3 * sizeof(constants.REAL_t))
point_new[0], point_new[1], point_new[2] = end
ray.push_back(point_new)
point = ray.back()
value = self.value(point)

while True:
gg = grad.value(point)
norm = sqrt(gg[0]**2 + gg[1]**2 + gg[2]**2)
for idx in range(3):
gg[idx] /= norm
if coord_sys == "spherical":
gg[1] /= point[0]
gg[2] /= point[0] * sin(point[1])
point_new = <constants.REAL_t *> malloc(3 * sizeof(constants.REAL_t))
for idx in range(3):
point_new[idx] = point[idx] - step_size * gg[idx]
point_1back = ray.back()
ray.push_back(point_new)
point = ray.back()
try:
value_1back = value
value = self.value(point)
if value_1back <= value or np.isnan(value):
break
except ValueError:
for idx in range(ray.size()-1):
free(ray[idx])
return (None)
ray_np = np.zeros((ray.size()-1, 3), dtype=constants.DTYPE_REAL)
for idx in range(ray.size()-1):
for jdx in range(3):
ray_np[idx, jdx] = ray[idx][jdx]
free(ray[idx])
return (np.flipud(ray_np))


cpdef constants.REAL_t value(
ScalarField3D self,
constants.REAL_t[:] point,
constants.REAL_t null=np.nan
):
"""
value(self, point, null=numpy.nan)
Expand Down
21 changes: 21 additions & 0 deletions pykonal/solver.pyx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,8 @@ instantiated as below:
solver = pykonal.EikonalSolver(coord_sys="cartesian")
"""

import warnings

# Python thrid-party imports.
import numpy as np

Expand Down Expand Up @@ -116,9 +118,16 @@ cdef class EikonalSolver(object):
@property
def norm(self):
"""
*DEPRECATED*
[*Read-only*, :class:`numpy.ndarray`\ (shape=(N0,N1,N2,3), dtype=numpy.float)] 4D array of scaling
factors for gradient operator.
"""

warning_message = "This attribute is deprecated and will go away in a "\
"future version of PyKonal. Use pykonal.solver.EikonalSolver"\
".traveltime.norm instead."
warnings.warn(warning_message, DeprecationWarning)

try:
return (np.asarray(self.cy_norm))
except AttributeError:
Expand All @@ -136,8 +145,15 @@ cdef class EikonalSolver(object):
@property
def step_size(self):
"""
*DEPRECATED*
[*Read only*, :class:`float`] Step size used for ray tracing.
"""

warning_message = "This attribute is deprecated and will go away in a "\
"future version of PyKonal. Use pykonal.solver.EikonalSolver"\
".traveltime.step_size instead."
warnings.warn(warning_message, DeprecationWarning)

return (self.norm[~np.isclose(self.norm, 0)].min() / 4)


Expand Down Expand Up @@ -411,6 +427,11 @@ cdef class EikonalSolver(object):
cdef fields.ScalarField3D traveltime
cdef str coord_sys

warning_message = "This method is deprecated and will disappear in a "\
"future version of PyKonal. Use pykonal.EikonalSolver.traveltime."\
"trace_ray() instead."
warnings.warn(warning_message, DeprecationWarning)


coord_sys = self.coord_sys
step_size = self.step_size
Expand Down

0 comments on commit e407e0c

Please sign in to comment.