[yt-svn] commit/yt: 12 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Thu Aug 20 09:17:23 PDT 2015
12 new commits in yt:
https://bitbucket.org/yt_analysis/yt/commits/ee9401b188eb/
Changeset: ee9401b188eb
Branch: yt
User: jisuoqing
Date: 2015-08-17 08:53:18+00:00
Summary: Add line integral convolution for plot modification callback
Affected #: 6 files
diff -r bc21601eb3e3d89810ca3b5e914797f934fe1ef8 -r ee9401b188eb934abceb245013796fbfd6db38ed doc/source/visualizing/callbacks.rst
--- a/doc/source/visualizing/callbacks.rst
+++ b/doc/source/visualizing/callbacks.rst
@@ -463,6 +463,34 @@
s.annotate_streamlines('velocity_x', 'velocity_y')
s.save()
+.. _annotate-line-integral-convolution:
+
+Overplot Line Integral Convolution
+~~~~~~~~~~~~~~~~~~~~
+
+.. function:: annotate_line_integral_convolution(self, field_x, field_y, \
+ texture=None, kernellen=50., \
+ lim=(0.5,0.6), cmap='binary', \
+ const_alpha=False, alpha=None)
+
+ (This is a proxy for
+ :class:`~yt.visualization.plot_modifications.LineIntegralConvolutionCallback`.)
+
+ Add line integral convolution to any plot, using the ``field_x`` and ``field_y``
+ from the associated data. A white noise background will be used for ``texture``
+ as default. Adjust ``lim`` to change the visibility of line integral convolution
+ plot. When ``const_alpha=False``, spatially varying alpha based on the values of
+ line integral convolution will be applied; otherwise a constant value of alpha
+ should be provided.
+
+.. python-script::
+
+ import yt
+ ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+ s = yt.SlicePlot(ds, 'z', 'density', center='c', width=(20, 'kpc'))
+ s.annotate_line_integral_convolution('velocity_x', 'velocity_y',lim=(0.5,0.65))
+ s.save()
+
.. _annotate-text:
Overplot Text
diff -r bc21601eb3e3d89810ca3b5e914797f934fe1ef8 -r ee9401b188eb934abceb245013796fbfd6db38ed yt/utilities/lib/api.py
--- a/yt/utilities/lib/api.py
+++ b/yt/utilities/lib/api.py
@@ -29,3 +29,4 @@
from .write_array import *
from .mesh_utilities import *
from .ContourFinding import *
+from .line_integral_convolution import *
diff -r bc21601eb3e3d89810ca3b5e914797f934fe1ef8 -r ee9401b188eb934abceb245013796fbfd6db38ed yt/utilities/lib/line_integral_convolution.pyx
--- /dev/null
+++ b/yt/utilities/lib/line_integral_convolution.pyx
@@ -0,0 +1,109 @@
+"""
+Utilities for line integral convolution annotation
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2015, yt Development Team.
+#
+# Code originally from Scipy Cookbook (http://wiki.scipy.org/Cookbook/LineIntegralConvolution),
+# with bug fixed which leads to crash when non equal-size vector field in two
+# dimensions is provided.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+
+import numpy as np
+cimport numpy as np
+
+cdef void _advance_2d(double vx, double vy,
+ int* x, int* y,
+ double* fx, double* fy,
+ int w, int h):
+ cdef double tx, ty
+ if vx>=0:
+ tx = (1-fx[0])/vx
+ else:
+ tx = -fx[0]/vx
+ if vy>=0:
+ ty = (1-fy[0])/vy
+ else:
+ ty = -fy[0]/vy
+ if tx<ty:
+ if vx>=0:
+ x[0]+=1
+ fx[0]=0
+ else:
+ x[0]-=1
+ fx[0]=1
+ fy[0]+=tx*vy
+ else:
+ if vy>=0:
+ y[0]+=1
+ fy[0]=0
+ else:
+ y[0]-=1
+ fy[0]=1
+ fx[0]+=ty*vx
+ if x[0]>=w:
+ x[0]=w-1
+ if x[0]<0:
+ x[0]=0
+ if y[0]<0:
+ y[0]=0
+ if y[0]>=h:
+ y[0]=h-1
+
+def line_integral_convolution_2d(
+ np.ndarray[double, ndim=3] vectors,
+ np.ndarray[double, ndim=2] texture,
+ np.ndarray[double, ndim=1] kernel):
+ cdef int i,j,l,x,y
+ cdef int h,w,kernellen
+ cdef int t
+ cdef double fx, fy, tx, ty
+ cdef np.ndarray[double, ndim=2] result
+
+ w = vectors.shape[0]
+ h = vectors.shape[1]
+ t = vectors.shape[2]
+
+ kernellen = kernel.shape[0]
+ result = np.zeros((w,h),dtype=np.double)
+
+ temp = np.copy(vectors[...,0])
+ vectors[...,0] = vectors[...,1]
+ vectors[...,1] = temp
+
+ for i in range(w):
+ for j in range(h):
+ x = i
+ y = j
+ fx = 0.5
+ fy = 0.5
+
+ l = kernellen//2
+ result[i,j] += kernel[l]*texture[x,y]
+ while l<kernellen-1:
+ _advance_2d(vectors[x,y,0],vectors[x,y,1],
+ &x, &y, &fx, &fy, w, h)
+ l+=1
+ result[i,j] += kernel[l]*texture[x,y]
+
+ x = i
+ y = j
+ fx = 0.5
+ fy = 0.5
+
+ while l>0:
+ _advance_2d(-vectors[x,y,0],-vectors[x,y,1],
+ &x, &y, &fx, &fy, w, h)
+ l-=1
+ result[i,j] += kernel[l]*texture[x,y]
+
+ return result
diff -r bc21601eb3e3d89810ca3b5e914797f934fe1ef8 -r ee9401b188eb934abceb245013796fbfd6db38ed yt/utilities/lib/setup.py
--- a/yt/utilities/lib/setup.py
+++ b/yt/utilities/lib/setup.py
@@ -155,6 +155,8 @@
config.add_extension("amr_kdtools",
["yt/utilities/lib/amr_kdtools.pyx"],
libraries=["m"], depends=["yt/utilities/lib/fp_utils.pxd"])
+ config.add_extension("line_integral_convolution",
+ ["yt/utilities/lib/line_integral_convolution.pyx"])
config.add_subpackage("tests")
if os.environ.get("GPERFTOOLS", "no").upper() != "NO":
diff -r bc21601eb3e3d89810ca3b5e914797f934fe1ef8 -r ee9401b188eb934abceb245013796fbfd6db38ed yt/visualization/plot_modifications.py
--- a/yt/visualization/plot_modifications.py
+++ b/yt/visualization/plot_modifications.py
@@ -23,6 +23,7 @@
from matplotlib.patches import Circle
from matplotlib.colors import colorConverter
+from matplotlib import cm
from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar
from yt.funcs import *
@@ -37,6 +38,8 @@
from yt.utilities.lib.geometry_utils import triangle_plane_intersect
from yt.analysis_modules.cosmological_observation.light_ray.light_ray \
import periodic_ray
+from yt.utilities.lib.line_integral_convolution \
+ import line_integral_convolution_2d
import warnings
from . import _MPL
@@ -2055,3 +2058,126 @@
lcb(plot)
return plot
+
+class LineIntegralConvolutionCallback(PlotCallback):
+ """
+ annotate_line_integral_convolution(field_x, field_y, texture=None,
+ kernellen=50., lim=(0.5,0.6),
+ cmap='binary', const_alpha=False,
+ alpha=None):
+
+ Add the line integral convolution to the plot for vector fields
+ visualization. Two component of vector fields needed to be provided
+ (i.e., velocity_x and velocity_y, magentic_field_x and magnetic_field_y).
+
+ Parameters
+ ----------
+
+ field_x, field_y : string
+ The names of two components of vector field which will be visualized
+
+ texture : 2-d array, optional
+ Texture will be convolved when computing line integral convolution.
+ A white noise background will be used as default.
+
+ kernellen : float, optional
+ The lens of kernel for convolution.
+
+ lim : 2-element tuple, list, or array, optional
+ The value of line integral convolution will be clipped to the range
+ of lim, which will enhance the visibility of the plot.
+
+ cmap : float, optional
+ The name of colormap for line integral convolution plot.
+
+ const_alpha : boolean, optional
+ If set to False (by default), spatially varying alpha based on the
+ values of line integral convolution will be applied; otherwise a
+ constant value of alpha is used.
+
+ alpha : float, optional
+ The alpha value for line integral convolution plot. Will be ignored
+ when const_alpha = False.
+
+ Example
+ -------
+
+ >>> import yt
+ >>> ds = yt.load('Enzo_64/DD0020/data0020')
+ >>> s = yt.SlicePlot(ds, 'z', 'density')
+ >>> s.annotate_line_integral_convolution('velocity_x','velocity_y')
+ """
+ _type_name = "line_integral_convolution"
+ def __init__(self, field_x, field_y, texture=None, kernellen=50.,
+ lim=(0.5,0.6), cmap='binary', const_alpha=False, alpha=None):
+ PlotCallback.__init__(self)
+ self.field_x = field_x
+ self.field_y = field_y
+ self.texture = texture
+ self.kernellen = kernellen
+ self.lim = lim
+ self.cmap = cmap
+ self.const_alpha = const_alpha
+ self.alpha = alpha
+
+ def __call__(self, plot):
+ x0, x1 = plot.xlim
+ y0, y1 = plot.ylim
+ xx0, xx1 = plot._axes.get_xlim()
+ yy0, yy1 = plot._axes.get_ylim()
+ extent = [xx0,xx1,yy0,yy1]
+
+ plot._axes.hold(True)
+ nx = plot.image._A.shape[0]
+ ny = plot.image._A.shape[1]
+ pixX = _MPL.Pixelize(plot.data['px'],
+ plot.data['py'],
+ plot.data['pdx'],
+ plot.data['pdy'],
+ plot.data[self.field_x],
+ int(nx), int(ny),
+ (x0, x1, y0, y1),).transpose()
+ pixY = _MPL.Pixelize(plot.data['px'],
+ plot.data['py'],
+ plot.data['pdx'],
+ plot.data['pdy'],
+ plot.data[self.field_y],
+ int(nx), int(ny),
+ (x0, x1, y0, y1),).transpose()
+ vectors = np.concatenate((pixX[...,np.newaxis],
+ pixY[...,np.newaxis]),axis=2)
+
+ if self.texture == None:
+ self.texture = np.random.rand(nx,ny).astype(np.double)
+ elif self.texture.shape != (nx,ny):
+ raise SyntaxError("'texture' must have the same shape "
+ "with that of output image (%d, %d)" % (nx,ny))
+
+ kernel = np.sin(np.arange(self.kernellen)*np.pi/self.kernellen)
+ kernel = kernel.astype(np.double)
+
+ lic_data = line_integral_convolution_2d(vectors,self.texture,kernel)
+ lic_data = np.flipud(lic_data / lic_data.max())
+ lic_data_clip = np.clip(lic_data,self.lim[0],self.lim[1])
+
+ if self.const_alpha:
+ if self.alpha == None:
+ self.alpha = 0.8
+ print("'const_alpha' is set to True but the value of alpha "
+ "is not provided; will use alpha = %f as default"
+ % self.alpha)
+ plot._axes.imshow(lic_data_clip, extent=extent, cmap=self.cmap,
+ alpha=self.alpha)
+ else:
+ if self.alpha != None:
+ print("The provided alpha value %f will be ignored when "
+ "'const_alpha' is set to False" % self.alpha)
+ lic_data_rgba = cm.ScalarMappable(norm=None, cmap=self.cmap).\
+ to_rgba(lic_data_clip)
+ lic_data_clip_rescale = (lic_data_clip - self.lim[0]) \
+ / (self.lim[1]-self.lim[0])
+ lic_data_rgba[...,3] = lic_data_clip_rescale
+ plot._axes.imshow(lic_data_rgba, extent=extent)
+ plot._axes.hold(False)
+
+ return plot
\ No newline at end of file
diff -r bc21601eb3e3d89810ca3b5e914797f934fe1ef8 -r ee9401b188eb934abceb245013796fbfd6db38ed yt/visualization/tests/test_callbacks.py
--- a/yt/visualization/tests/test_callbacks.py
+++ b/yt/visualization/tests/test_callbacks.py
@@ -52,6 +52,7 @@
# X scale
# material_boundary
# X ray
+# X
@contextlib.contextmanager
def _cleanup_fname():
@@ -326,3 +327,25 @@
max_level=3, cmap="gist_stern")
p.save(prefix)
+def test_line_integral_convolution_callback():
+ with _cleanup_fname() as prefix:
+ ds = fake_amr_ds(fields =
+ ("density", "velocity_x", "velocity_y", "velocity_z"))
+ for ax in 'xyz':
+ p = ProjectionPlot(ds, ax, "density")
+ p.annotate_line_integral_convolution("velocity_x", "velocity_y")
+ yield assert_fname, p.save(prefix)[0]
+ p = ProjectionPlot(ds, ax, "density", weight_field="density")
+ p.annotate_line_integral_convolution("velocity_x", "velocity_y")
+ yield assert_fname, p.save(prefix)[0]
+ p = SlicePlot(ds, ax, "density")
+ p.annotate_line_integral_convolution("velocity_x", "velocity_y")
+ yield assert_fname, p.save(prefix)[0]
+ # Now we'll check a few additional minor things
+ p = SlicePlot(ds, "x", "density")
+ p.annotate_line_integral_convolution("velocity_x", "velocity_y",
+ kernellen=100., lim=(0.4,0.7),
+ cmap='algae', const_alpha=True,
+ alpha=0.9)
+ p.save(prefix)
+
https://bitbucket.org/yt_analysis/yt/commits/f1e72eef2dc5/
Changeset: f1e72eef2dc5
Branch: yt
User: jisuoqing
Date: 2015-08-17 08:55:52+00:00
Summary: A minor fix
Affected #: 1 file
diff -r ee9401b188eb934abceb245013796fbfd6db38ed -r f1e72eef2dc5f16021c5e43451755a227c41d83a yt/visualization/tests/test_callbacks.py
--- a/yt/visualization/tests/test_callbacks.py
+++ b/yt/visualization/tests/test_callbacks.py
@@ -52,7 +52,7 @@
# X scale
# material_boundary
# X ray
-# X
+# X line_integral_convolution
@contextlib.contextmanager
def _cleanup_fname():
https://bitbucket.org/yt_analysis/yt/commits/e13cac97d439/
Changeset: e13cac97d439
Branch: yt
User: jisuoqing
Date: 2015-08-17 09:13:17+00:00
Summary: Minor change
Affected #: 2 files
diff -r f1e72eef2dc5f16021c5e43451755a227c41d83a -r e13cac97d43965661cc29cc845738b879ac30d1c doc/source/visualizing/callbacks.rst
--- a/doc/source/visualizing/callbacks.rst
+++ b/doc/source/visualizing/callbacks.rst
@@ -488,7 +488,7 @@
import yt
ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
s = yt.SlicePlot(ds, 'z', 'density', center='c', width=(20, 'kpc'))
- s.annotate_line_integral_convolution('velocity_x', 'velocity_y',lim=(0.5,0.65))
+ s.annotate_line_integral_convolution('velocity_x', 'velocity_y', lim=(0.5,0.65))
s.save()
.. _annotate-text:
diff -r f1e72eef2dc5f16021c5e43451755a227c41d83a -r e13cac97d43965661cc29cc845738b879ac30d1c yt/visualization/plot_modifications.py
--- a/yt/visualization/plot_modifications.py
+++ b/yt/visualization/plot_modifications.py
@@ -2105,7 +2105,8 @@
>>> import yt
>>> ds = yt.load('Enzo_64/DD0020/data0020')
>>> s = yt.SlicePlot(ds, 'z', 'density')
- >>> s.annotate_line_integral_convolution('velocity_x','velocity_y')
+ >>> s.annotate_line_integral_convolution('velocity_x', 'velocity_y',\
+ lim=(0.5,0.65))
"""
_type_name = "line_integral_convolution"
def __init__(self, field_x, field_y, texture=None, kernellen=50.,
@@ -2180,4 +2181,4 @@
plot._axes.imshow(lic_data_rgba, extent=extent)
plot._axes.hold(False)
- return plot
\ No newline at end of file
+ return plot
https://bitbucket.org/yt_analysis/yt/commits/cdeea919652c/
Changeset: cdeea919652c
Branch: yt
User: jisuoqing
Date: 2015-08-17 09:39:53+00:00
Summary: Enhance the flexibility of varying alpha option
Affected #: 2 files
diff -r e13cac97d43965661cc29cc845738b879ac30d1c -r cdeea919652cb4f5b7ec955187f28f8f0e7a5449 doc/source/visualizing/callbacks.rst
--- a/doc/source/visualizing/callbacks.rst
+++ b/doc/source/visualizing/callbacks.rst
@@ -471,17 +471,18 @@
.. function:: annotate_line_integral_convolution(self, field_x, field_y, \
texture=None, kernellen=50., \
lim=(0.5,0.6), cmap='binary', \
- const_alpha=False, alpha=None)
+ alpha=0.8, const_alpha=False)
(This is a proxy for
:class:`~yt.visualization.plot_modifications.LineIntegralConvolutionCallback`.)
Add line integral convolution to any plot, using the ``field_x`` and ``field_y``
from the associated data. A white noise background will be used for ``texture``
- as default. Adjust ``lim`` to change the visibility of line integral convolution
- plot. When ``const_alpha=False``, spatially varying alpha based on the values of
- line integral convolution will be applied; otherwise a constant value of alpha
- should be provided.
+ as default. Adjust the bounds of ``lim`` in the range of ``[0, 1]`` which applies
+ upper and lower bounds to the values of line integral convolution and enhance
+ the visibility of plots. When ``const_alpha=False``, alpha will be weighted
+ spatially by the values of line integral convolution; otherwise a constant value
+ of the given alpha is used.
.. python-script::
diff -r e13cac97d43965661cc29cc845738b879ac30d1c -r cdeea919652cb4f5b7ec955187f28f8f0e7a5449 yt/visualization/plot_modifications.py
--- a/yt/visualization/plot_modifications.py
+++ b/yt/visualization/plot_modifications.py
@@ -2061,10 +2061,10 @@
class LineIntegralConvolutionCallback(PlotCallback):
"""
- annotate_line_integral_convolution(field_x, field_y, texture=None,
- kernellen=50., lim=(0.5,0.6),
- cmap='binary', const_alpha=False,
- alpha=None):
+ annotate_line_integral_convolution(field_x, field_y, texture=None,
+ kernellen=50., lim=(0.5,0.6),
+ cmap='binary', alpha=0.8,
+ const_alpha=False):
Add the line integral convolution to the plot for vector fields
visualization. Two component of vector fields needed to be provided
@@ -2084,20 +2084,21 @@
The lens of kernel for convolution.
lim : 2-element tuple, list, or array, optional
- The value of line integral convolution will be clipped to the range
- of lim, which will enhance the visibility of the plot.
+ The value of line integral convolution will be clipped to the range
+ of lim, which applies upper and lower bounds to the values of line
+ integral convolution and enhance the visibility of plots. Each element
+ should be in the range of [0,1].
cmap : float, optional
The name of colormap for line integral convolution plot.
+ alpha : float, optional
+ The alpha value for line integral convolution plot.
+
const_alpha : boolean, optional
- If set to False (by default), spatially varying alpha based on the
- values of line integral convolution will be applied; otherwise a
- constant value of alpha is used.
-
- alpha : float, optional
- The alpha value for line integral convolution plot. Will be ignored
- when const_alpha = False.
+ If set to False (by default), alpha will be weighted spatially by
+ the values of line integral convolution; otherwise a constant value
+ of the given alpha is used.
Example
-------
@@ -2110,7 +2111,7 @@
"""
_type_name = "line_integral_convolution"
def __init__(self, field_x, field_y, texture=None, kernellen=50.,
- lim=(0.5,0.6), cmap='binary', const_alpha=False, alpha=None):
+ lim=(0.5,0.6), cmap='binary', alpha=0.8, const_alpha=False):
PlotCallback.__init__(self)
self.field_x = field_x
self.field_y = field_y
@@ -2118,8 +2119,8 @@
self.kernellen = kernellen
self.lim = lim
self.cmap = cmap
+ self.alpha = alpha
self.const_alpha = const_alpha
- self.alpha = alpha
def __call__(self, plot):
x0, x1 = plot.xlim
@@ -2162,22 +2163,14 @@
lic_data_clip = np.clip(lic_data,self.lim[0],self.lim[1])
if self.const_alpha:
- if self.alpha == None:
- self.alpha = 0.8
- print("'const_alpha' is set to True but the value of alpha "
- "is not provided; will use alpha = %f as default"
- % self.alpha)
plot._axes.imshow(lic_data_clip, extent=extent, cmap=self.cmap,
alpha=self.alpha)
else:
- if self.alpha != None:
- print("The provided alpha value %f will be ignored when "
- "'const_alpha' is set to False" % self.alpha)
lic_data_rgba = cm.ScalarMappable(norm=None, cmap=self.cmap).\
to_rgba(lic_data_clip)
lic_data_clip_rescale = (lic_data_clip - self.lim[0]) \
/ (self.lim[1]-self.lim[0])
- lic_data_rgba[...,3] = lic_data_clip_rescale
+ lic_data_rgba[...,3] = lic_data_clip_rescale * self.alpha
plot._axes.imshow(lic_data_rgba, extent=extent)
plot._axes.hold(False)
https://bitbucket.org/yt_analysis/yt/commits/4aa017ec1039/
Changeset: 4aa017ec1039
Branch: yt
User: jisuoqing
Date: 2015-08-17 09:46:10+00:00
Summary: Minor change in test script
Affected #: 1 file
diff -r cdeea919652cb4f5b7ec955187f28f8f0e7a5449 -r 4aa017ec10390ac2fcb865b4910b4cc90281d4a8 yt/visualization/tests/test_callbacks.py
--- a/yt/visualization/tests/test_callbacks.py
+++ b/yt/visualization/tests/test_callbacks.py
@@ -345,7 +345,7 @@
p = SlicePlot(ds, "x", "density")
p.annotate_line_integral_convolution("velocity_x", "velocity_y",
kernellen=100., lim=(0.4,0.7),
- cmap='algae', const_alpha=True,
- alpha=0.9)
+ cmap='algae', alpha=0.9,
+ const_alpha=True)
p.save(prefix)
https://bitbucket.org/yt_analysis/yt/commits/066585a9e095/
Changeset: 066585a9e095
Branch: yt
User: jisuoqing
Date: 2015-08-17 18:19:25+00:00
Summary: Minor fix on Docs
Affected #: 1 file
diff -r 4aa017ec10390ac2fcb865b4910b4cc90281d4a8 -r 066585a9e095de8c7924f1c975140f3d8e7ff69f doc/source/visualizing/callbacks.rst
--- a/doc/source/visualizing/callbacks.rst
+++ b/doc/source/visualizing/callbacks.rst
@@ -466,7 +466,7 @@
.. _annotate-line-integral-convolution:
Overplot Line Integral Convolution
-~~~~~~~~~~~~~~~~~~~~
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
.. function:: annotate_line_integral_convolution(self, field_x, field_y, \
texture=None, kernellen=50., \
https://bitbucket.org/yt_analysis/yt/commits/67235bf40559/
Changeset: 67235bf40559
Branch: yt
User: jisuoqing
Date: 2015-08-17 18:51:37+00:00
Summary: Typo fixed
Affected #: 1 file
diff -r 066585a9e095de8c7924f1c975140f3d8e7ff69f -r 67235bf40559cf2b38f11f73b57b15b03074c02b yt/visualization/plot_modifications.py
--- a/yt/visualization/plot_modifications.py
+++ b/yt/visualization/plot_modifications.py
@@ -2089,7 +2089,7 @@
integral convolution and enhance the visibility of plots. Each element
should be in the range of [0,1].
- cmap : float, optional
+ cmap : string, optional
The name of colormap for line integral convolution plot.
alpha : float, optional
https://bitbucket.org/yt_analysis/yt/commits/6df0176fd795/
Changeset: 6df0176fd795
Branch: yt
User: jisuoqing
Date: 2015-08-17 20:40:03+00:00
Summary: Switched to ds.coordinates.pixelize
Affected #: 1 file
diff -r 67235bf40559cf2b38f11f73b57b15b03074c02b -r 6df0176fd79552e3fed92344f01c64e394b382c2 yt/visualization/plot_modifications.py
--- a/yt/visualization/plot_modifications.py
+++ b/yt/visualization/plot_modifications.py
@@ -2127,25 +2127,23 @@
y0, y1 = plot.ylim
xx0, xx1 = plot._axes.get_xlim()
yy0, yy1 = plot._axes.get_ylim()
+ bounds = [x0,x1,y0,y1]
extent = [xx0,xx1,yy0,yy1]
plot._axes.hold(True)
nx = plot.image._A.shape[0]
ny = plot.image._A.shape[1]
- pixX = _MPL.Pixelize(plot.data['px'],
- plot.data['py'],
- plot.data['pdx'],
- plot.data['pdy'],
- plot.data[self.field_x],
- int(nx), int(ny),
- (x0, x1, y0, y1),).transpose()
- pixY = _MPL.Pixelize(plot.data['px'],
- plot.data['py'],
- plot.data['pdx'],
- plot.data['pdy'],
- plot.data[self.field_y],
- int(nx), int(ny),
- (x0, x1, y0, y1),).transpose()
+ pixX = plot.data.ds.coordinates.pixelize(plot.data.axis,
+ plot.data,
+ self.field_x,
+ bounds,
+ (nx,ny))
+ pixY = plot.data.ds.coordinates.pixelize(plot.data.axis,
+ plot.data,
+ self.field_y,
+ bounds,
+ (nx,ny))
+
vectors = np.concatenate((pixX[...,np.newaxis],
pixY[...,np.newaxis]),axis=2)
https://bitbucket.org/yt_analysis/yt/commits/206e32b3c6bb/
Changeset: 206e32b3c6bb
Branch: yt
User: jisuoqing
Date: 2015-08-17 23:07:26+00:00
Summary: Update docstring
Affected #: 1 file
diff -r 6df0176fd79552e3fed92344f01c64e394b382c2 -r 206e32b3c6bb4193124ae96995f42460474db95d yt/visualization/plot_modifications.py
--- a/yt/visualization/plot_modifications.py
+++ b/yt/visualization/plot_modifications.py
@@ -2076,12 +2076,14 @@
field_x, field_y : string
The names of two components of vector field which will be visualized
- texture : 2-d array, optional
+ texture : 2-d array with the same shape of image, optional
Texture will be convolved when computing line integral convolution.
A white noise background will be used as default.
kernellen : float, optional
- The lens of kernel for convolution.
+ The lens of kernel for convolution, which is the length over which the
+ convolution will be performed. For longer kernellen, longer streamline
+ structure will appear.
lim : 2-element tuple, list, or array, optional
The value of line integral convolution will be clipped to the range
@@ -2167,7 +2169,7 @@
lic_data_rgba = cm.ScalarMappable(norm=None, cmap=self.cmap).\
to_rgba(lic_data_clip)
lic_data_clip_rescale = (lic_data_clip - self.lim[0]) \
- / (self.lim[1]-self.lim[0])
+ / (self.lim[1] - self.lim[0])
lic_data_rgba[...,3] = lic_data_clip_rescale * self.alpha
plot._axes.imshow(lic_data_rgba, extent=extent)
plot._axes.hold(False)
https://bitbucket.org/yt_analysis/yt/commits/4daf3abec54c/
Changeset: 4daf3abec54c
Branch: yt
User: jisuoqing
Date: 2015-08-19 18:20:08+00:00
Summary: Minor changes
Affected #: 2 files
diff -r 206e32b3c6bb4193124ae96995f42460474db95d -r 4daf3abec54c1bb12f415741c238367b56ffa690 yt/utilities/lib/line_integral_convolution.pyx
--- a/yt/utilities/lib/line_integral_convolution.pyx
+++ b/yt/utilities/lib/line_integral_convolution.pyx
@@ -20,7 +20,9 @@
import numpy as np
cimport numpy as np
+cimport cython
+ at cython.cdivision(True)
cdef void _advance_2d(double vx, double vy,
int* x, int* y,
double* fx, double* fy,
@@ -76,9 +78,7 @@
kernellen = kernel.shape[0]
result = np.zeros((w,h),dtype=np.double)
- temp = np.copy(vectors[...,0])
- vectors[...,0] = vectors[...,1]
- vectors[...,1] = temp
+ vectors = vectors[...,::-1].copy()
for i in range(w):
for j in range(h):
diff -r 206e32b3c6bb4193124ae96995f42460474db95d -r 4daf3abec54c1bb12f415741c238367b56ffa690 yt/visualization/plot_modifications.py
--- a/yt/visualization/plot_modifications.py
+++ b/yt/visualization/plot_modifications.py
@@ -2167,7 +2167,7 @@
alpha=self.alpha)
else:
lic_data_rgba = cm.ScalarMappable(norm=None, cmap=self.cmap).\
- to_rgba(lic_data_clip)
+ to_rgba(lic_data_clip)
lic_data_clip_rescale = (lic_data_clip - self.lim[0]) \
/ (self.lim[1] - self.lim[0])
lic_data_rgba[...,3] = lic_data_clip_rescale * self.alpha
https://bitbucket.org/yt_analysis/yt/commits/d2e9e1f113a4/
Changeset: d2e9e1f113a4
Branch: yt
User: jisuoqing
Date: 2015-08-20 00:13:50+00:00
Summary: Merged yt_analysis/yt into yt
Affected #: 12 files
diff -r 4daf3abec54c1bb12f415741c238367b56ffa690 -r d2e9e1f113a46b53c031f109eb2003fc0f748d4f doc/source/examining/loading_data.rst
--- a/doc/source/examining/loading_data.rst
+++ b/doc/source/examining/loading_data.rst
@@ -771,7 +771,7 @@
.. code-block:: python
- from yt.frontends.sph.definitions import gadget_field_specs
+ from yt.frontends.gadget.definitions import gadget_field_specs
gadget_field_specs["my_field_def"] = my_field_def
Please also feel free to issue a pull request with any new field
@@ -871,7 +871,7 @@
----------------
See :ref:`loading-numpy-array` and
-:func:`~yt.frontends.sph.data_structures.load_amr_grids` for more detail.
+:func:`~yt.frontends.stream.data_structures.load_amr_grids` for more detail.
It is possible to create native yt dataset from Python's dictionary
that describes set of rectangular patches of data of possibly varying
diff -r 4daf3abec54c1bb12f415741c238367b56ffa690 -r d2e9e1f113a46b53c031f109eb2003fc0f748d4f yt/data_objects/construction_data_containers.py
--- a/yt/data_objects/construction_data_containers.py
+++ b/yt/data_objects/construction_data_containers.py
@@ -682,11 +682,13 @@
def RightEdge(self):
return self.right_edge
- def deposit(self, positions, fields = None, method = None):
+ def deposit(self, positions, fields = None, method = None,
+ kernel_name = 'cubic'):
cls = getattr(particle_deposit, "deposit_%s" % method, None)
if cls is None:
raise YTParticleDepositionNotImplemented(method)
- op = cls(self.ActiveDimensions.prod()) # We allocate number of zones, not number of octs
+ # We allocate number of zones, not number of octs
+ op = cls(self.ActiveDimensions.prod(), kernel_name)
op.initialize()
op.process_grid(self, positions, fields)
vals = op.finalize()
@@ -1787,5 +1789,3 @@
else:
mylog.error("Problem uploading.")
return upload_id
-
-
diff -r 4daf3abec54c1bb12f415741c238367b56ffa690 -r d2e9e1f113a46b53c031f109eb2003fc0f748d4f yt/data_objects/grid_patch.py
--- a/yt/data_objects/grid_patch.py
+++ b/yt/data_objects/grid_patch.py
@@ -330,12 +330,14 @@
def particle_operation(self, *args, **kwargs):
raise NotImplementedError
- def deposit(self, positions, fields = None, method = None):
+ def deposit(self, positions, fields = None, method = None,
+ kernel_name = 'cubic'):
# Here we perform our particle deposition.
cls = getattr(particle_deposit, "deposit_%s" % method, None)
if cls is None:
raise YTParticleDepositionNotImplemented(method)
- op = cls(self.ActiveDimensions.prod()) # We allocate number of zones, not number of octs
+ # We allocate number of zones, not number of octs
+ op = cls(self.ActiveDimensions.prod(), kernel_name)
op.initialize()
op.process_grid(self, positions, fields)
vals = op.finalize()
diff -r 4daf3abec54c1bb12f415741c238367b56ffa690 -r d2e9e1f113a46b53c031f109eb2003fc0f748d4f yt/data_objects/octree_subset.py
--- a/yt/data_objects/octree_subset.py
+++ b/yt/data_objects/octree_subset.py
@@ -141,7 +141,8 @@
self._domain_ind = di
return self._domain_ind
- def deposit(self, positions, fields = None, method = None):
+ def deposit(self, positions, fields = None, method = None,
+ kernel_name='cubic'):
r"""Operate on the mesh, in a particle-against-mesh fashion, with
exclusively local input.
@@ -176,7 +177,8 @@
raise YTParticleDepositionNotImplemented(method)
nz = self.nz
nvals = (nz, nz, nz, (self.domain_ind >= 0).sum())
- op = cls(nvals) # We allocate number of zones, not number of octs
+ # We allocate number of zones, not number of octs
+ op = cls(nvals, kernel_name)
op.initialize()
mylog.debug("Depositing %s (%s^3) particles into %s Octs",
positions.shape[0], positions.shape[0]**0.3333333, nvals[-1])
@@ -192,7 +194,8 @@
return np.asfortranarray(vals)
def smooth(self, positions, fields = None, index_fields = None,
- method = None, create_octree = False, nneighbors = 64):
+ method = None, create_octree = False, nneighbors = 64,
+ kernel_name = 'cubic'):
r"""Operate on the mesh, in a particle-against-mesh fashion, with
non-local input.
@@ -258,7 +261,7 @@
nz = self.nz
mdom_ind = self.domain_ind
nvals = (nz, nz, nz, (mdom_ind >= 0).sum())
- op = cls(nvals, len(fields), nneighbors)
+ op = cls(nvals, len(fields), nneighbors, kernel_name)
op.initialize()
mylog.debug("Smoothing %s particles into %s Octs",
positions.shape[0], nvals[-1])
@@ -280,7 +283,7 @@
return vals
def particle_operation(self, positions, fields = None,
- method = None, nneighbors = 64):
+ method = None, nneighbors = 64, kernel_name = 'cubic'):
r"""Operate on particles, in a particle-against-particle fashion.
This uses the octree indexing system to call a "smoothing" operation
@@ -335,7 +338,7 @@
nz = self.nz
mdom_ind = self.domain_ind
nvals = (nz, nz, nz, (mdom_ind >= 0).sum())
- op = cls(nvals, len(fields), nneighbors)
+ op = cls(nvals, len(fields), nneighbors, kernel_name)
op.initialize()
mylog.debug("Smoothing %s particles into %s Octs",
positions.shape[0], nvals[-1])
diff -r 4daf3abec54c1bb12f415741c238367b56ffa690 -r d2e9e1f113a46b53c031f109eb2003fc0f748d4f yt/data_objects/unstructured_mesh.py
--- a/yt/data_objects/unstructured_mesh.py
+++ b/yt/data_objects/unstructured_mesh.py
@@ -105,13 +105,15 @@
def select_tcoords(self, dobj):
raise NotImplementedError
- def deposit(self, positions, fields = None, method = None):
+ def deposit(self, positions, fields = None, method = None,
+ kernel_name = 'cubic'):
raise NotImplementedError
# Here we perform our particle deposition.
cls = getattr(particle_deposit, "deposit_%s" % method, None)
if cls is None:
raise YTParticleDepositionNotImplemented(method)
- op = cls(self.ActiveDimensions.prod()) # We allocate number of zones, not number of octs
+ # We allocate number of zones, not number of octs
+ op = cls(self.ActiveDimensions.prod(), kernel_name)
op.initialize()
op.process_grid(self, positions, fields)
vals = op.finalize()
@@ -200,4 +202,3 @@
else:
self._last_count = mask.sum()
return mask
-
diff -r 4daf3abec54c1bb12f415741c238367b56ffa690 -r d2e9e1f113a46b53c031f109eb2003fc0f748d4f yt/fields/particle_fields.py
--- a/yt/fields/particle_fields.py
+++ b/yt/fields/particle_fields.py
@@ -766,8 +766,12 @@
def add_volume_weighted_smoothed_field(ptype, coord_name, mass_name,
smoothing_length_name, density_name, smoothed_field, registry,
- nneighbors = None):
- field_name = ("deposit", "%s_smoothed_%s" % (ptype, smoothed_field))
+ nneighbors = None, kernel_name = 'cubic'):
+ if kernel_name == 'cubic':
+ field_name = ("deposit", "%s_smoothed_%s" % (ptype, smoothed_field))
+ else:
+ field_name = ("deposit", "%s_%s_smoothed_%s" % (ptype, kernel_name,
+ smoothed_field))
field_units = registry[ptype, smoothed_field].units
def _vol_weight(field, data):
pos = data[ptype, coord_name]
@@ -789,7 +793,8 @@
# volume_weighted smooth operations return lists of length 1.
rv = data.smooth(pos, [mass, hsml, dens, quan],
method="volume_weighted",
- create_octree=True)[0]
+ create_octree=True,
+ kernel_name=kernel_name)[0]
rv[np.isnan(rv)] = 0.0
# Now some quick unit conversions.
rv = data.apply_units(rv, field_units)
@@ -816,8 +821,12 @@
units = "code_length")
return [field_name]
-def add_density_kernel(ptype, coord_name, mass_name, registry, nneighbors = 64):
- field_name = (ptype, "smoothed_density")
+def add_density_kernel(ptype, coord_name, mass_name, registry, nneighbors = 64,
+ kernel_name = 'cubic'):
+ if kernel_name == 'cubic':
+ field_name = (ptype, "smoothed_density")
+ else:
+ field_name = (ptype, "%s_smoothed_density" % (kernel_name))
field_units = registry[ptype, mass_name].units
def _nth_neighbor(field, data):
pos = data[ptype, coord_name]
@@ -827,7 +836,8 @@
densities = mass * 0.0
data.particle_operation(pos, [mass, densities],
method="density",
- nneighbors = nneighbors)
+ nneighbors = nneighbors,
+ kernel_name = kernel_name)
ones = pos.prod(axis=1) # Get us in code_length**3
ones[:] = 1.0
densities /= ones
diff -r 4daf3abec54c1bb12f415741c238367b56ffa690 -r d2e9e1f113a46b53c031f109eb2003fc0f748d4f yt/frontends/artio/data_structures.py
--- a/yt/frontends/artio/data_structures.py
+++ b/yt/frontends/artio/data_structures.py
@@ -133,7 +133,8 @@
tr = dict((field, v) for field, v in zip(fields, tr))
return tr
- def deposit(self, positions, fields = None, method = None):
+ def deposit(self, positions, fields = None, method = None,
+ kernel_name = 'cubic'):
# Here we perform our particle deposition.
if fields is None: fields = []
cls = getattr(particle_deposit, "deposit_%s" % method, None)
@@ -141,7 +142,8 @@
raise YTParticleDepositionNotImplemented(method)
nz = self.nz
nvals = (nz, nz, nz, self.ires.size)
- op = cls(nvals) # We allocate number of zones, not number of octs
+ # We allocate number of zones, not number of octs
+ op = cls(nvals, kernel_name)
op.initialize()
mylog.debug("Depositing %s (%s^3) particles into %s Root Mesh",
positions.shape[0], positions.shape[0]**0.3333333, nvals[-1])
diff -r 4daf3abec54c1bb12f415741c238367b56ffa690 -r d2e9e1f113a46b53c031f109eb2003fc0f748d4f yt/frontends/artio/setup.py
--- a/yt/frontends/artio/setup.py
+++ b/yt/frontends/artio/setup.py
@@ -19,7 +19,8 @@
depends=artio_sources +
["yt/utilities/lib/fp_utils.pxd",
"yt/geometry/oct_container.pxd",
- "yt/geometry/selection_routines.pxd"])
+ "yt/geometry/selection_routines.pxd",
+ "yt/geometry/particle_deposit.pxd"])
config.make_config_py() # installs __config__.py
#config.make_svn_version_py()
return config
diff -r 4daf3abec54c1bb12f415741c238367b56ffa690 -r d2e9e1f113a46b53c031f109eb2003fc0f748d4f yt/geometry/particle_deposit.pxd
--- a/yt/geometry/particle_deposit.pxd
+++ b/yt/geometry/particle_deposit.pxd
@@ -38,7 +38,7 @@
# Standard SPH kernel for use with the Grid method #
####################################################
-cdef inline np.float64_t sph_kernel(np.float64_t x) nogil:
+cdef inline np.float64_t sph_kernel_cubic(np.float64_t x) nogil:
cdef np.float64_t kernel
if x <= 0.5:
kernel = 1.-6.*x*x*(1.-x)
@@ -48,8 +48,55 @@
kernel = 0.
return kernel
+########################################################
+# Alternative SPH kernels for use with the Grid method #
+########################################################
+
+# quartic spline
+cdef inline np.float64_t sph_kernel_quartic(np.float64_t x):
+ cdef np.float64_t kernel
+ cdef np.float64_t C = 5.**6/512/np.pi
+ if x < 1:
+ kernel = (1.-x)**4
+ if x < 3./5:
+ kernel -= 5*(3./5-x)**4
+ if x < 1./5:
+ kernel += 10*(1./5-x)**4
+ else:
+ kernel = 0.
+ return kernel * C
+
+# quintic spline
+cdef inline np.float64_t sph_kernel_quintic(np.float64_t x):
+ cdef np.float64_t kernel
+ cdef np.float64_t C = 3.**7/40/np.pi
+ if x < 1:
+ kernel = (1.-x)**5
+ if x < 2./3:
+ kernel -= 6*(2./3-x)**5
+ if x < 1./3:
+ kernel += 15*(1./3-x)**5
+ else:
+ kernel = 0.
+ return kernel * C
+
+# I don't know the way to use a dict in a cdef class.
+# So in order to mimic a registry functionality,
+# I manually created a function to lookup the kernel functions.
+ctypedef np.float64_t (*kernel_func) (np.float64_t)
+cdef inline kernel_func get_kernel_func(str kernel_name):
+ if kernel_name == 'cubic':
+ return sph_kernel_cubic
+ elif kernel_name == 'quartic':
+ return sph_kernel_quartic
+ elif kernel_name == 'quintic':
+ return sph_kernel_quintic
+ else:
+ raise NotImplementedError
+
cdef class ParticleDepositOperation:
# We assume each will allocate and define their own temporary storage
+ cdef kernel_func sph_kernel
cdef public object nvals
cdef public int update_values
cdef void process(self, int dim[3], np.float64_t left_edge[3],
diff -r 4daf3abec54c1bb12f415741c238367b56ffa690 -r d2e9e1f113a46b53c031f109eb2003fc0f748d4f yt/geometry/particle_deposit.pyx
--- a/yt/geometry/particle_deposit.pyx
+++ b/yt/geometry/particle_deposit.pyx
@@ -25,9 +25,10 @@
OctreeContainer, OctInfo
cdef class ParticleDepositOperation:
- def __init__(self, nvals):
+ def __init__(self, nvals, kernel_name):
self.nvals = nvals
self.update_values = 0 # This is the default
+ self.sph_kernel = get_kernel_func(kernel_name)
def initialize(self, *args):
raise NotImplementedError
@@ -227,7 +228,7 @@
dist = idist[0] + idist[1] + idist[2]
# Calculate distance in multiples of the smoothing length
dist = sqrt(dist) / fields[0]
- self.temp[gind(i,j,k,dim) + offset] = sph_kernel(dist)
+ self.temp[gind(i,j,k,dim) + offset] = self.sph_kernel(dist)
kernel_sum += self.temp[gind(i,j,k,dim) + offset]
# Having found the kernel, deposit accordingly into gdata
for i from ib0[0] <= i <= ib1[0]:
@@ -493,4 +494,3 @@
return self.onnfield
deposit_nearest = NNParticleField
-
diff -r 4daf3abec54c1bb12f415741c238367b56ffa690 -r d2e9e1f113a46b53c031f109eb2003fc0f748d4f yt/geometry/particle_smooth.pxd
--- a/yt/geometry/particle_smooth.pxd
+++ b/yt/geometry/particle_smooth.pxd
@@ -22,7 +22,7 @@
from fp_utils cimport *
from oct_container cimport Oct, OctAllocationContainer, OctreeContainer
-from .particle_deposit cimport sph_kernel, gind
+from .particle_deposit cimport kernel_func, get_kernel_func, gind
cdef extern from "platform_dep.h":
void *alloca(int)
@@ -34,6 +34,7 @@
cdef class ParticleSmoothOperation:
# We assume each will allocate and define their own temporary storage
+ cdef kernel_func sph_kernel
cdef public object nvals
cdef np.float64_t DW[3]
cdef int nfields
diff -r 4daf3abec54c1bb12f415741c238367b56ffa690 -r d2e9e1f113a46b53c031f109eb2003fc0f748d4f yt/geometry/particle_smooth.pyx
--- a/yt/geometry/particle_smooth.pyx
+++ b/yt/geometry/particle_smooth.pyx
@@ -74,7 +74,7 @@
opos[2] = ipos[2]
cdef class ParticleSmoothOperation:
- def __init__(self, nvals, nfields, max_neighbors):
+ def __init__(self, nvals, nfields, max_neighbors, kernel_name):
# This is the set of cells, in grids, blocks or octs, we are handling.
cdef int i
self.nvals = nvals
@@ -83,6 +83,7 @@
self.neighbors = <NeighborList *> malloc(
sizeof(NeighborList) * self.maxn)
self.neighbor_reset()
+ self.sph_kernel = get_kernel_func(kernel_name)
def initialize(self, *args):
raise NotImplementedError
@@ -630,7 +631,7 @@
# Usually this density has been computed
dens = fields[2][pn]
if dens == 0.0: continue
- weight = mass * sph_kernel(sqrt(r2) / hsml) / dens
+ weight = mass * self.sph_kernel(sqrt(r2) / hsml) / dens
# Mass of the particle times the value
for fi in range(self.nfields - 3):
val = fields[fi + 3][pn]
@@ -756,7 +757,7 @@
for pn in range(self.curn):
mass = fields[0][self.neighbors[pn].pn]
r2 = self.neighbors[pn].r2
- lw = sph_kernel(sqrt(r2) / hsml)
+ lw = self.sph_kernel(sqrt(r2) / hsml)
dens += mass * lw
weight = (4.0/3.0) * 3.1415926 * hsml**3
fields[1][offset] = dens/weight
https://bitbucket.org/yt_analysis/yt/commits/98b9f2c2d791/
Changeset: 98b9f2c2d791
Branch: yt
User: MatthewTurk
Date: 2015-08-20 16:17:12+00:00
Summary: Merged in jisuoqing/yt (pull request #1702)
[NEW] [DOCS] Add line integral convolution for plot modification callback
Affected #: 6 files
diff -r 1e76952e015dfb14db5d97333c182edc658d4aaa -r 98b9f2c2d791d921eccc8e2edc419d863eb14650 doc/source/visualizing/callbacks.rst
--- a/doc/source/visualizing/callbacks.rst
+++ b/doc/source/visualizing/callbacks.rst
@@ -463,6 +463,35 @@
s.annotate_streamlines('velocity_x', 'velocity_y')
s.save()
+.. _annotate-line-integral-convolution:
+
+Overplot Line Integral Convolution
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. function:: annotate_line_integral_convolution(self, field_x, field_y, \
+ texture=None, kernellen=50., \
+ lim=(0.5,0.6), cmap='binary', \
+ alpha=0.8, const_alpha=False)
+
+ (This is a proxy for
+ :class:`~yt.visualization.plot_modifications.LineIntegralConvolutionCallback`.)
+
+ Add line integral convolution to any plot, using the ``field_x`` and ``field_y``
+ from the associated data. A white noise background will be used for ``texture``
+ as default. Adjust the bounds of ``lim`` in the range of ``[0, 1]`` which applies
+ upper and lower bounds to the values of line integral convolution and enhance
+ the visibility of plots. When ``const_alpha=False``, alpha will be weighted
+ spatially by the values of line integral convolution; otherwise a constant value
+ of the given alpha is used.
+
+.. python-script::
+
+ import yt
+ ds = yt.load("IsolatedGalaxy/galaxy0030/galaxy0030")
+ s = yt.SlicePlot(ds, 'z', 'density', center='c', width=(20, 'kpc'))
+ s.annotate_line_integral_convolution('velocity_x', 'velocity_y', lim=(0.5,0.65))
+ s.save()
+
.. _annotate-text:
Overplot Text
diff -r 1e76952e015dfb14db5d97333c182edc658d4aaa -r 98b9f2c2d791d921eccc8e2edc419d863eb14650 yt/utilities/lib/api.py
--- a/yt/utilities/lib/api.py
+++ b/yt/utilities/lib/api.py
@@ -29,3 +29,4 @@
from .write_array import *
from .mesh_utilities import *
from .ContourFinding import *
+from .line_integral_convolution import *
diff -r 1e76952e015dfb14db5d97333c182edc658d4aaa -r 98b9f2c2d791d921eccc8e2edc419d863eb14650 yt/utilities/lib/line_integral_convolution.pyx
--- /dev/null
+++ b/yt/utilities/lib/line_integral_convolution.pyx
@@ -0,0 +1,109 @@
+"""
+Utilities for line integral convolution annotation
+
+
+
+"""
+
+#-----------------------------------------------------------------------------
+# Copyright (c) 2015, yt Development Team.
+#
+# Code originally from Scipy Cookbook (http://wiki.scipy.org/Cookbook/LineIntegralConvolution),
+# with bug fixed which leads to crash when non equal-size vector field in two
+# dimensions is provided.
+#
+# Distributed under the terms of the Modified BSD License.
+#
+# The full license is in the file COPYING.txt, distributed with this software.
+#-----------------------------------------------------------------------------
+
+
+import numpy as np
+cimport numpy as np
+cimport cython
+
+ at cython.cdivision(True)
+cdef void _advance_2d(double vx, double vy,
+ int* x, int* y,
+ double* fx, double* fy,
+ int w, int h):
+ cdef double tx, ty
+ if vx>=0:
+ tx = (1-fx[0])/vx
+ else:
+ tx = -fx[0]/vx
+ if vy>=0:
+ ty = (1-fy[0])/vy
+ else:
+ ty = -fy[0]/vy
+ if tx<ty:
+ if vx>=0:
+ x[0]+=1
+ fx[0]=0
+ else:
+ x[0]-=1
+ fx[0]=1
+ fy[0]+=tx*vy
+ else:
+ if vy>=0:
+ y[0]+=1
+ fy[0]=0
+ else:
+ y[0]-=1
+ fy[0]=1
+ fx[0]+=ty*vx
+ if x[0]>=w:
+ x[0]=w-1
+ if x[0]<0:
+ x[0]=0
+ if y[0]<0:
+ y[0]=0
+ if y[0]>=h:
+ y[0]=h-1
+
+def line_integral_convolution_2d(
+ np.ndarray[double, ndim=3] vectors,
+ np.ndarray[double, ndim=2] texture,
+ np.ndarray[double, ndim=1] kernel):
+ cdef int i,j,l,x,y
+ cdef int h,w,kernellen
+ cdef int t
+ cdef double fx, fy, tx, ty
+ cdef np.ndarray[double, ndim=2] result
+
+ w = vectors.shape[0]
+ h = vectors.shape[1]
+ t = vectors.shape[2]
+
+ kernellen = kernel.shape[0]
+ result = np.zeros((w,h),dtype=np.double)
+
+ vectors = vectors[...,::-1].copy()
+
+ for i in range(w):
+ for j in range(h):
+ x = i
+ y = j
+ fx = 0.5
+ fy = 0.5
+
+ l = kernellen//2
+ result[i,j] += kernel[l]*texture[x,y]
+ while l<kernellen-1:
+ _advance_2d(vectors[x,y,0],vectors[x,y,1],
+ &x, &y, &fx, &fy, w, h)
+ l+=1
+ result[i,j] += kernel[l]*texture[x,y]
+
+ x = i
+ y = j
+ fx = 0.5
+ fy = 0.5
+
+ while l>0:
+ _advance_2d(-vectors[x,y,0],-vectors[x,y,1],
+ &x, &y, &fx, &fy, w, h)
+ l-=1
+ result[i,j] += kernel[l]*texture[x,y]
+
+ return result
diff -r 1e76952e015dfb14db5d97333c182edc658d4aaa -r 98b9f2c2d791d921eccc8e2edc419d863eb14650 yt/utilities/lib/setup.py
--- a/yt/utilities/lib/setup.py
+++ b/yt/utilities/lib/setup.py
@@ -155,6 +155,8 @@
config.add_extension("amr_kdtools",
["yt/utilities/lib/amr_kdtools.pyx"],
libraries=["m"], depends=["yt/utilities/lib/fp_utils.pxd"])
+ config.add_extension("line_integral_convolution",
+ ["yt/utilities/lib/line_integral_convolution.pyx"])
config.add_subpackage("tests")
if os.environ.get("GPERFTOOLS", "no").upper() != "NO":
diff -r 1e76952e015dfb14db5d97333c182edc658d4aaa -r 98b9f2c2d791d921eccc8e2edc419d863eb14650 yt/visualization/plot_modifications.py
--- a/yt/visualization/plot_modifications.py
+++ b/yt/visualization/plot_modifications.py
@@ -23,6 +23,7 @@
from matplotlib.patches import Circle
from matplotlib.colors import colorConverter
+from matplotlib import cm
from mpl_toolkits.axes_grid1.anchored_artists import AnchoredSizeBar
from yt.funcs import *
@@ -37,6 +38,8 @@
from yt.utilities.lib.geometry_utils import triangle_plane_intersect
from yt.analysis_modules.cosmological_observation.light_ray.light_ray \
import periodic_ray
+from yt.utilities.lib.line_integral_convolution \
+ import line_integral_convolution_2d
import warnings
from . import _MPL
@@ -2055,3 +2058,120 @@
lcb(plot)
return plot
+
+class LineIntegralConvolutionCallback(PlotCallback):
+ """
+ annotate_line_integral_convolution(field_x, field_y, texture=None,
+ kernellen=50., lim=(0.5,0.6),
+ cmap='binary', alpha=0.8,
+ const_alpha=False):
+
+ Add the line integral convolution to the plot for vector fields
+ visualization. Two component of vector fields needed to be provided
+ (i.e., velocity_x and velocity_y, magentic_field_x and magnetic_field_y).
+
+ Parameters
+ ----------
+
+ field_x, field_y : string
+ The names of two components of vector field which will be visualized
+
+ texture : 2-d array with the same shape of image, optional
+ Texture will be convolved when computing line integral convolution.
+ A white noise background will be used as default.
+
+ kernellen : float, optional
+ The lens of kernel for convolution, which is the length over which the
+ convolution will be performed. For longer kernellen, longer streamline
+ structure will appear.
+
+ lim : 2-element tuple, list, or array, optional
+ The value of line integral convolution will be clipped to the range
+ of lim, which applies upper and lower bounds to the values of line
+ integral convolution and enhance the visibility of plots. Each element
+ should be in the range of [0,1].
+
+ cmap : string, optional
+ The name of colormap for line integral convolution plot.
+
+ alpha : float, optional
+ The alpha value for line integral convolution plot.
+
+ const_alpha : boolean, optional
+ If set to False (by default), alpha will be weighted spatially by
+ the values of line integral convolution; otherwise a constant value
+ of the given alpha is used.
+
+ Example
+ -------
+
+ >>> import yt
+ >>> ds = yt.load('Enzo_64/DD0020/data0020')
+ >>> s = yt.SlicePlot(ds, 'z', 'density')
+ >>> s.annotate_line_integral_convolution('velocity_x', 'velocity_y',\
+ lim=(0.5,0.65))
+ """
+ _type_name = "line_integral_convolution"
+ def __init__(self, field_x, field_y, texture=None, kernellen=50.,
+ lim=(0.5,0.6), cmap='binary', alpha=0.8, const_alpha=False):
+ PlotCallback.__init__(self)
+ self.field_x = field_x
+ self.field_y = field_y
+ self.texture = texture
+ self.kernellen = kernellen
+ self.lim = lim
+ self.cmap = cmap
+ self.alpha = alpha
+ self.const_alpha = const_alpha
+
+ def __call__(self, plot):
+ x0, x1 = plot.xlim
+ y0, y1 = plot.ylim
+ xx0, xx1 = plot._axes.get_xlim()
+ yy0, yy1 = plot._axes.get_ylim()
+ bounds = [x0,x1,y0,y1]
+ extent = [xx0,xx1,yy0,yy1]
+
+ plot._axes.hold(True)
+ nx = plot.image._A.shape[0]
+ ny = plot.image._A.shape[1]
+ pixX = plot.data.ds.coordinates.pixelize(plot.data.axis,
+ plot.data,
+ self.field_x,
+ bounds,
+ (nx,ny))
+ pixY = plot.data.ds.coordinates.pixelize(plot.data.axis,
+ plot.data,
+ self.field_y,
+ bounds,
+ (nx,ny))
+
+ vectors = np.concatenate((pixX[...,np.newaxis],
+ pixY[...,np.newaxis]),axis=2)
+
+ if self.texture == None:
+ self.texture = np.random.rand(nx,ny).astype(np.double)
+ elif self.texture.shape != (nx,ny):
+ raise SyntaxError("'texture' must have the same shape "
+ "with that of output image (%d, %d)" % (nx,ny))
+
+ kernel = np.sin(np.arange(self.kernellen)*np.pi/self.kernellen)
+ kernel = kernel.astype(np.double)
+
+ lic_data = line_integral_convolution_2d(vectors,self.texture,kernel)
+ lic_data = np.flipud(lic_data / lic_data.max())
+ lic_data_clip = np.clip(lic_data,self.lim[0],self.lim[1])
+
+ if self.const_alpha:
+ plot._axes.imshow(lic_data_clip, extent=extent, cmap=self.cmap,
+ alpha=self.alpha)
+ else:
+ lic_data_rgba = cm.ScalarMappable(norm=None, cmap=self.cmap).\
+ to_rgba(lic_data_clip)
+ lic_data_clip_rescale = (lic_data_clip - self.lim[0]) \
+ / (self.lim[1] - self.lim[0])
+ lic_data_rgba[...,3] = lic_data_clip_rescale * self.alpha
+ plot._axes.imshow(lic_data_rgba, extent=extent)
+ plot._axes.hold(False)
+
+ return plot
diff -r 1e76952e015dfb14db5d97333c182edc658d4aaa -r 98b9f2c2d791d921eccc8e2edc419d863eb14650 yt/visualization/tests/test_callbacks.py
--- a/yt/visualization/tests/test_callbacks.py
+++ b/yt/visualization/tests/test_callbacks.py
@@ -52,6 +52,7 @@
# X scale
# material_boundary
# X ray
+# X line_integral_convolution
@contextlib.contextmanager
def _cleanup_fname():
@@ -326,3 +327,25 @@
max_level=3, cmap="gist_stern")
p.save(prefix)
+def test_line_integral_convolution_callback():
+ with _cleanup_fname() as prefix:
+ ds = fake_amr_ds(fields =
+ ("density", "velocity_x", "velocity_y", "velocity_z"))
+ for ax in 'xyz':
+ p = ProjectionPlot(ds, ax, "density")
+ p.annotate_line_integral_convolution("velocity_x", "velocity_y")
+ yield assert_fname, p.save(prefix)[0]
+ p = ProjectionPlot(ds, ax, "density", weight_field="density")
+ p.annotate_line_integral_convolution("velocity_x", "velocity_y")
+ yield assert_fname, p.save(prefix)[0]
+ p = SlicePlot(ds, ax, "density")
+ p.annotate_line_integral_convolution("velocity_x", "velocity_y")
+ yield assert_fname, p.save(prefix)[0]
+ # Now we'll check a few additional minor things
+ p = SlicePlot(ds, "x", "density")
+ p.annotate_line_integral_convolution("velocity_x", "velocity_y",
+ kernellen=100., lim=(0.4,0.7),
+ cmap='algae', alpha=0.9,
+ const_alpha=True)
+ p.save(prefix)
+
Repository URL: https://bitbucket.org/yt_analysis/yt/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the yt-svn
mailing list