ENH/DEP: use more a more robust, external implementation of the Line Integral Convolution (LIC) algorithm

Author: neutrinoceros

.gitignore

@@ -54,7 +54,6 @@ yt/utilities/lib/image_utilities.c
 yt/utilities/lib/interpolators.c
 yt/utilities/lib/kdtree.c
 yt/utilities/lib/lenses.c
-yt/utilities/lib/line_integral_convolution.c
 yt/utilities/lib/embree_mesh/mesh_construction.cpp
 yt/utilities/lib/embree_mesh/mesh_intersection.cpp
 yt/utilities/lib/embree_mesh/mesh_samplers.cpp

minimal_requirements.txt

@@ -47,6 +47,8 @@ pyparsing==3.1.4
     #   packaging
 python-dateutil==2.9.0.post0
     # via matplotlib
+rlic==0.3.2
+    # via yt (pyproject.toml)
 six==1.16.0
     # via python-dateutil
 sympy==1.13.2

pyproject.toml

@@ -51,6 +51,7 @@ dependencies = [
     "numpy!=2.0.1 ; platform_machine=='arm64' and platform_system=='Darwin'",
     "packaging>=20.9",
     "pillow>=8.3.2",
+    "rlic>=0.3.2",
     "tomli-w>=0.4.0",
     "tqdm>=3.4.0",
     "unyt>=2.9.2",
@@ -432,7 +433,6 @@ ignore = [
   "yt/utilities/lib/image_utilities.c",
   "yt/utilities/lib/interpolators.c",
   "yt/utilities/lib/lenses.c",
-  "yt/utilities/lib/line_integral_convolution.c",
   "yt/utilities/lib/marching_cubes.cpp",
   "yt/utilities/lib/mesh_triangulation.c",
   "yt/utilities/lib/mesh_utilities.c",

yt/utilities/lib/api.py

@@ -5,7 +5,6 @@
 from .grid_traversal import *
 from .image_utilities import *
 from .interpolators import *
-from .line_integral_convolution import *
 from .marching_cubes import *
 from .mesh_utilities import *
 from .misc_utilities import *

yt/utilities/lib/line_integral_convolution.pyx

@@ -9,6 +9,7 @@
 
 import matplotlib
 import numpy as np
+import rlic
 from unyt import unyt_quantity
 
 from yt._maintenance.deprecation import issue_deprecation_warning
@@ -30,7 +31,6 @@
     YTUnsupportedPlotCallback,
 )
 from yt.utilities.lib.geometry_utils import triangle_plane_intersect
-from yt.utilities.lib.line_integral_convolution import line_integral_convolution_2d
 from yt.utilities.lib.mesh_triangulation import triangulate_indices
 from yt.utilities.lib.pixelization_routines import (
     pixelize_cartesian,
@@ -3238,10 +3238,10 @@ def __init__(
         field_x,
         field_y,
         texture=None,
-        kernellen=50.0,
+        kernellen=50,
         lim=(0.5, 0.6),
         cmap="binary",
-        alpha=0.8,
+        alpha=0.5,
         const_alpha=False,
     ):
         self.field_x = field_x
@@ -3270,29 +3270,18 @@ def __call__(self, plot):
             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 is None:
-            prng = np.random.RandomState(0x4D3D3D3)
-            self.texture = prng.random_sample((nx, ny))
-        elif self.texture.shape != (nx, ny):
-            raise ValueError(
-                "'texture' must have the same shape "
-                f"with that of output image ({nx}, {ny})"
-            )
+            prng = np.random.default_rng(0x4D3D3D3)
+            self.texture = prng.random((nx, ny))
 
         kernel = np.sin(
             np.arange(self.kernellen, dtype="float64") * np.pi / self.kernellen
         )
 
-        lic_data = line_integral_convolution_2d(vectors, self.texture, kernel)
+        lic_data = rlic.convolve(self.texture, pixX, pixY, kernel=kernel)
         lic_data = lic_data / lic_data.max()
         lic_data_clip = np.clip(lic_data, self.lim[0], self.lim[1])
 
-        mask = ~(np.isfinite(pixX) & np.isfinite(pixY))
-        lic_data[mask] = np.nan
-        lic_data_clip[mask] = np.nan
-
         if plot._swap_axes:
             lic_data_clip = lic_data_clip.transpose()
             extent = (extent[2], extent[3], extent[0], extent[1])