Added additional volume rendering

docs/source/visualisation/volume_render.rst

@@ -155,6 +155,122 @@ the above example is shown below.
    # A 256 x 256 x 256 cube with dimensions of temperature
    temp_cube = mass_weighted_temp_cube / mass_cube
 
+
+Rendering
+---------
+
+We provide a volume rendering function that can be used to make images highlighting
+specific density contours. The notable function here is
+:meth:``swiftsimio.visualisation.volume_render.visualise_render``. This takes
+in your volume rendering, along with a colour map and centers, to create
+these highlights. The example below shows how to use this.
+
+.. code-block:: python
+
+   import matplotlib.pyplot as plt
+   import numpy as np
+   from matplotlib.colors import LogNorm
+   
+   from swiftsimio import load
+   from swiftsimio.visualisation import volume_render
+   
+   # Load the data
+   data = load("test_data/eagle_6.hdf5")
+   
+   # Rough location of an interesting galaxy in the volume.
+   region = [
+       0.225 * data.metadata.boxsize[0],
+       0.275 * data.metadata.boxsize[0],
+       0.12 * data.metadata.boxsize[1],
+       0.17 * data.metadata.boxsize[1],
+       0.45 * data.metadata.boxsize[2],
+       0.5 * data.metadata.boxsize[2],
+   ]
+   
+   # Render the volume (note 1024 is reasonably high resolution so this won't complete
+   # immediately; you should consider using 256, etc. for testing).
+   rendered = volume_render.render_gas(data, resolution=1024, region=region, parallel=True)
+   
+   # Quick view! By projecting along the final axis you can get
+   # the projected density from the rendered image.
+   plt.imsave("volume_render_quick_view.png", LogNorm()(rendered.sum(-1)))
+
+Here we can see the quick view of this image. It's just a regular density projection:
+
+.. image:: volume_render_quick_view.png
+
+.. code-block:: python
+   
+   # Now we will move onto the real volume rendering. Let's use the log of the density;
+   # using the real density leads to low contrast images.
+   log_rendered = np.log10(rendered)
+   
+   # The volume rendering function expects centers of 'bins' and widths. These
+   # bins actually represent gaussian functions around a specific density (or other
+   # visualization quantity). The brightest pixel value is at center. We will
+   # visualise this later!
+   width = 0.1
+   std = np.std(log_rendered)
+   mean = np.mean(log_rendered)
+   
+   # It's helpful to choose the centers relative to the data you have. When making
+   # a movie, you will obviously want to choose the centers to be the same for each
+   # frame.
+   centers = [mean + x * std for x in [1.0, 3.0, 5.0, 7.0]]
+   
+   # This will visualize your render options. The centers are shown as gaussians and
+   # vertical lines.
+   fig, ax = volume_render.visualise_render_options(
+       centers=centers, widths=width, cmap="viridis"
+   )
+   
+   histogram, edges = np.histogram(
+       log_rendered.flat,
+       bins=128,
+       range=(min(centers) - 5.0 * width, max(centers) + 5.0 * width),
+   )
+   bc = (edges[:-1] + edges[1:]) / 2.0
+   
+   # The normalization here is the height of a gaussian!
+   ax.plot(bc, histogram / (np.max(histogram) * np.sqrt(2.0 * np.pi) * width))
+   ax.semilogy()
+   ax.set_xlabel("$\\log_{10}(\\rho)$")
+   
+   plt.savefig("volume_render_options.png")
+
+This function :meth:`swiftsimio.visualisation.volume_render.visualise_render_options` allows
+you to see what densities your rendering is picking out:
+
+.. image:: volume_render_options.png
+
+.. code-block:: python   
+   # Now we can really visualize the rendering.
+   img, norms = volume_render.visualise_render(
+       log_rendered,
+       centers,
+       widths=width,
+       cmap="viridis",
+   )
+   
+   # Sometimes, these images can be a bit dark. You can increase the brightness using
+   # tools like PIL or in your favourite image editor.
+   from PIL import Image, ImageEnhance
+   
+   pilimg = Image.fromarray((img * 255.0).astype(np.uint8))
+   enhanced = ImageEnhance.Contrast(ImageEnhance.Brightness(pilimg).enhance(2.0)).enhance(
+       1.2
+   )
+   
+   enhanced.save("volume_render_example.png")
+
+Which produces the image:
+
+.. image:: volume_render_example.png
+
+Once you have this base image, you can always use your photo editor to tweak it further.
+In particular, open the 'levels' panel and play around with the sliders!
+
+
 Lower-level API
 ---------------
 

swiftsimio/optional_packages.py

@@ -39,6 +39,15 @@ def tqdm(x, *args, **kwargs):
     ASTROPY_AVAILABLE = False
 
 
+# matplotlib
+try:
+    import matplotlib.pyplot as plt
+
+    MATPLOTLIB_AVAILABLE = True
+except (ImportError, ModuleNotFoundError):
+    plt = None
+    MATPLOTLIB_AVAILABLE = False
+
 # Numba/CUDA
 try:
     from numba.cuda.cudadrv.error import CudaSupportError
@@ -83,3 +92,5 @@ def x(func):
 
     # For additional CUDA API access
     cuda = None
+
+

swiftsimio/visualisation/ray_trace.py

@@ -16,7 +16,7 @@
 from swiftsimio.reader import __SWIFTParticleDataset, SWIFTDataset
 from swiftsimio.visualisation.projection_backends.kernels import kernel_gamma, kernel_double_precision as kernel
 
-from swiftsimio.accelerated import jit
+from swiftsimio.accelerated import jit, prange, NUM_THREADS
 import unyt
 
 @jit(nopython=True, fastmath=True)
@@ -163,6 +163,52 @@ def core_panels(
 
     return output
 
+@jit(nopython=True, fastmath=True)
+def core_panels_parallel(
+    x: np.float64,
+    y: np.float64,
+    z: np.float64,
+    h: np.float32,
+    m: np.float32,
+    res: int,
+    panels: int,
+    min_z: np.float64,
+    max_z: np.float64,
+):
+    # Same as scatter, but executes in parallel! This is actually trivial,
+    # we just make NUM_THREADS images and add them together at the end.
+
+    number_of_particles = x.size
+    core_particles = number_of_particles // NUM_THREADS
+
+    output = np.zeros((res, res, panels), dtype=np.float32)
+
+    for thread in prange(NUM_THREADS):
+        # Left edge is easy, just start at 0 and go to 'final'
+        left_edge = thread * core_particles
+
+        # Right edge is harder in case of left over particles...
+        right_edge = thread + 1
+
+        if right_edge == NUM_THREADS:
+            right_edge = number_of_particles
+        else:
+            right_edge *= core_particles
+
+        output += core_panels(
+            x[left_edge:right_edge],
+            y[left_edge:right_edge],
+            z[left_edge:right_edge],
+            h[left_edge:right_edge],
+            m[left_edge:right_edge],
+            res,
+            panels,
+            min_z,
+            max_z,
+        )
+
+    return output
+
 
 def panel_pixel_grid(
     data: __SWIFTParticleDataset,