Merge pull request #150 from SWIFTSIM/complete-projections-new

Made sure visualisation functions include periodic copies of particles (attempt 2)

docs/source/visualisation/index.rst

@@ -17,6 +17,7 @@ common interface:
        project="masses", # Variable to project, e.g. masses, temperatures, etc.
        parallel=False, # Construct the image in (thread) parallel?
        region=None, # None, or a list telling which region to render_func_name
+       periodic=True, # Whether or not to apply periodic boundary conditions
    )
 
 The output of these functions comes with associated units and has the correct

docs/source/visualisation/projection.rst

@@ -30,7 +30,13 @@ Example
 
    # This creates a grid that has units msun / Mpc^2, and can be transformed like
    # any other unyt quantity
-   mass_map = project_gas(data, resolution=1024, project="masses", parallel=True)
+   mass_map = project_gas(
+       data,
+       resolution=1024,
+       project="masses",
+       parallel=True,
+       periodic=True,
+   )
 
    # Let's say we wish to save it as msun / kpc^2,
    from unyt import msun, kpc
@@ -61,13 +67,20 @@ this:
    data.gas.mass_weighted_temps = data.gas.masses * data.gas.temperatures
 
    # Map in msun / mpc^2
-   mass_map = project_gas(data, resolution=1024, project="masses", parallel=True)
+   mass_map = project_gas(
+       data,
+       resolution=1024,
+       project="masses",
+       parallel=True,
+       periodic=True,
+   )
    # Map in msun * K / mpc^2
    mass_weighted_temp_map = project_gas(
        data,
        resolution=1024,
        project="mass_weighted_temps",
-       parallel=True
+       parallel=True,
+       periodic=True,
    )
 
    temp_map = mass_weighted_temp_map / mass_map
@@ -132,13 +145,38 @@ Example:
       resolution=1024,
       project="entropies",
       parallel=True,
-      backend="subsampled"
+      backend="subsampled",
+      periodic=True,
    )
 
 This will likely look very similar to the image that you make with the default
 ``backend="fast"``, but will have a well-converged distribution at any resolution
 level.
 
+Periodic boundaries
+-------------------
+
+Cosmological simulations and many other simulations use periodic boundary
+conditions. This has implications for the particles at the edge of the
+simulation box: they can contribute to pixels on multiple sides of the image.
+If this effect is not taken into account, then the pixels close to the edge
+will have values that are too low because of missing contributions.
+
+All visualisation functions by default assume a periodic box. Rather than
+simply projecting each individual particle once, four additional periodic copies
+of each particle are also projected. Most copies will project outside the valid
+pixel range, but the copies that do not ensure that pixels close to the edge
+receive all necessary contributions. Thanks to Numba optimisations, the overhead
+of these additional copies is relatively small.
+
+There are some caveats with this approach. If you try to visualise a subset of
+the particles in the box (e.g. using a mask), then only periodic copies of
+particles in this subset will be used. If the subset does not include particles
+on the other side of the periodic boundary, then these will still be missing
+from the projection. The same is true if you visualise a region of the box.
+The periodic boundary wrapping is also not compatible with rotations (see below)
+and should therefore not be used together with a rotation.
+
 Rotations
 ---------
 
@@ -214,7 +252,13 @@ is shown in the ``velociraptor`` section.
 
    # Use project_gas_pixel_grid to generate projected images
 
-   common_arguments = dict(data=data, resolution=512, parallel=True, region=visualise_region)
+   common_arguments = dict(
+       data=data,
+       resolution=512,
+       parallel=True,
+       region=visualise_region,
+       periodic=False, # disable periodic boundaries when using rotations
+   )
 
    un_rotated = project_gas_pixel_grid(**common_arguments)
 
@@ -280,7 +324,8 @@ mass density map for dark matter. We provide a utility to do this through
        resolution=1024,
        project="masses",
        parallel=True,
-       region=None
+       region=None,
+       periodic=True,
    )
 
    from matplotlib.pyplot import imsave
@@ -322,6 +367,9 @@ To use this function, you will need:
 + Smoothing lengths for all particles, ``h``.
 + The resolution you wish to make your square image at, ``res``.
 
+Optionally, you will also need:
++ the size of the simulation box in x and y, ``box_x`` and ``box_y``.
+
 The key here is that only particles in the domain [0, 1] in x, and [0, 1] in y
 will be visible in the image. You may have particles outside of this range;
 they will not crash the code, and may even contribute to the image if their
@@ -338,3 +386,9 @@ such that it lives within this range. Then you may use the function as follows:
 ``out`` will be a 2D :mod:`numpy` grid of shape ``[res, res]``. You will need
 to re-scale this back to your original dimensions to get it in the correct units,
 and do not forget that it now represents the smoothed quantity per surface area.
+
+If the optional arguments ``box_x`` and ``box_y`` are provided, they should
+contain the simulation box size in the same re-scaled coordinates as ``x`` and
+``y``. The projection backend will then correctly apply periodic boundary
+wrapping. If ``box_x`` and ``box_y`` are not provided or set to 0, no
+periodic boundaries are applied.

docs/source/visualisation/slice.rst

@@ -36,7 +36,8 @@ Example
        z_slice=0.5 * data.metadata.boxsize[2],
        resolution=1024,
        project="masses",
-       parallel=True
+       parallel=True,
+       periodic=True,
    )
 
    # Let's say we wish to save it as g / cm^2,
@@ -72,7 +73,8 @@ in units of K / kpc^3 and we just want K) by dividing out by this:
        z_slice=0.5 * data.metadata.boxsize[2],
        resolution=1024,
        project="masses",
-       parallel=True
+       parallel=True,
+       periodic=True,
    )
 
    # Map in msun * K / mpc^3
@@ -81,7 +83,8 @@ in units of K / kpc^3 and we just want K) by dividing out by this:
        z_slice=0.5 * data.metadata.boxsize[2],
        resolution=1024,
        project="mass_weighted_temps",
-       parallel=True
+       parallel=True,
+       periodic=True,
    )
 
    temp_map = mass_weighted_temp_map / mass_map
@@ -101,6 +104,31 @@ loading data section should look something like:
 
 .. image:: temp_slice.png
 
+Periodic boundaries
+-------------------
+
+Cosmological simulations and many other simulations use periodic boundary
+conditions. This has implications for the particles at the edge of the
+simulation box: they can contribute to pixels on multiple sides of the image.
+If this effect is not taken into account, then the pixels close to the edge
+will have values that are too low because of missing contributions.
+
+All visualisation functions by default assume a periodic box. Rather than
+simply summing each individual particle once, eight additional periodic copies
+of each particle are also accounted for. Most copies will contribute outside the
+valid pixel range, but the copies that do not ensure that pixels close to the
+edge receive all necessary contributions. Thanks to Numba optimisations, the
+overhead of these additional copies is relatively small.
+
+There are some caveats with this approach. If you try to visualise a subset of
+the particles in the box (e.g. using a mask), then only periodic copies of
+particles in this subset will be used. If the subset does not include particles
+on the other side of the periodic boundary, then these will still be missing
+from the slice. The same is true if you visualise a region of the box.
+The periodic boundary wrapping is also not compatible with rotations (see below)
+and should therefore not be used together with a rotation.
+
+
 Rotations
 ---------
 
@@ -141,7 +169,8 @@ the above example is shown below.
        project="masses",
        rotation_matrix=matrix,
        rotation_center=center,
-       parallel=True
+       parallel=True,
+       periodic=False, # disable periodic boundaries when using rotations
    )
    
    # Map in msun * K / mpc^3
@@ -152,7 +181,8 @@ the above example is shown below.
        project="mass_weighted_temps",
        rotation_matrix=matrix,
        rotation_center=center,
-       parallel=True
+       parallel=True,
+       periodic=False,
    )
 
    temp_map = mass_weighted_temp_map / mass_map
@@ -192,6 +222,9 @@ To use this function, you will need:
 + Smoothing lengths for all particles, ``h``.
 + The resolution you wish to make your square image at, ``res``.
 
+Optionally, you will also need:
++ the size of the simulation box in x, y and z, ``box_x``, ``box_y`` and ``box_z``.
+
 The key here is that only particles in the domain [0, 1] in x and y will be
 visible in the image. You may have particles outside of this range; they will
 not crash the code, and may even contribute to the image if their smoothing
@@ -210,3 +243,9 @@ not need to lie in the domain [0, 1]. Then you may use the function as follows:
 ``out`` will be a 2D :mod:`numpy` grid of shape ``[res, res]``. You will need
 to re-scale this back to your original dimensions to get it in the correct units,
 and do not forget that it now represents the smoothed quantity per volume.
+
+If the optional arguments ``box_x``, ``box_y`` and ``box_z`` are provided, they
+should contain the simulation box size in the same re-scaled coordinates as 
+``x``, ``y`` and ``z``. The slicing function will then correctly apply
+periodic boundary wrapping. If ``box_x``, ``box_y`` and ``box_z`` are not
+provided or set to 0, no periodic boundaries are applied.

docs/source/visualisation/volume_render.rst

@@ -34,7 +34,8 @@ Example
        data,
        resolution=256,
        project="masses",
-       parallel=True
+       parallel=True,
+       periodic=True,
    )
 
 This basic demonstration creates a mass density cube.
@@ -59,20 +60,46 @@ this:
        data,
        resolution=256,
        project="masses",
-       parallel=True
+       parallel=True,
+       periodic=True,
    )
 
    # Map in msun * K / mpc^3
    mass_weighted_temp_cube = render_gas(
        data,
        resolution=256,
        project="mass_weighted_temps",
-       parallel=True
+       parallel=True,
+       periodic=True,
    )
 
    # A 256 x 256 x 256 cube with dimensions of temperature
    temp_cube = mass_weighted_temp_cube / mass_cube
 
+Periodic boundaries
+-------------------
+
+Cosmological simulations and many other simulations use periodic boundary
+conditions. This has implications for the particles at the edge of the
+simulation box: they can contribute to voxels on multiple sides of the image.
+If this effect is not taken into account, then the voxels close to the edge
+will have values that are too low because of missing contributions.
+
+All visualisation functions by default assume a periodic box. Rather than
+simply summing each individual particle once, eight additional periodic copies
+of each particle are also taken into account. Most copies will contribute
+outside the valid voxel range, but the copies that do not ensure that voxels
+close to the edge receive all necessary contributions. Thanks to Numba
+optimisations, the overhead of these additional copies is relatively small.
+
+There are some caveats with this approach. If you try to visualise a subset of
+the particles in the box (e.g. using a mask), then only periodic copies of
+particles in this subset will be used. If the subset does not include particles
+on the other side of the periodic boundary, then these will still be missing
+from the voxel cube. The same is true if you visualise a region of the box.
+The periodic boundary wrapping is also not compatible with rotations (see below)
+and should therefore not be used together with a rotation.
+
 Rotations
 ---------
 
@@ -110,7 +137,8 @@ the above example is shown below.
        project="masses",
        rotation_matrix=matrix,
        rotation_center=center,
-       parallel=True
+       parallel=True,
+       periodic=False, # disable periodic boundaries for rotations
    )
    
    # Map in msun * K / mpc^3
@@ -120,7 +148,8 @@ the above example is shown below.
        project="mass_weighted_temps",
        rotation_matrix=matrix,
        rotation_center=center,
-       parallel=True
+       parallel=True,
+       periodic=False,
    )
 
    # A 256 x 256 x 256 cube with dimensions of temperature
@@ -150,6 +179,9 @@ To use this function, you will need:
 + Smoothing lengths for all particles, ``h``.
 + The resolution you wish to make your cube at, ``res``.
 
+Optionally, you will also need:
++ the size of the simulation box in x, y and z, ``box_x``, ``box_y`` and ``box_z``.
+
 The key here is that only particles in the domain [0, 1] in x, [0, 1] in y,
 and [0, 1] in z. will be visible in the cube. You may have particles outside
 of this range; they will not crash the code, and may even contribute to the
@@ -168,3 +200,9 @@ function as follows:
 need to re-scale this back to your original dimensions to get it in the
 correct units, and do not forget that it now represents the smoothed quantity
 per volume.
+
+If the optional arguments ``box_x``, ``box_y`` and ``box_z`` are provided, they
+should contain the simulation box size in the same re-scaled coordinates as 
+``x``, ``y`` and ``z``. The rendering function will then correctly apply
+periodic boundary wrapping. If ``box_x``, ``box_y`` and ``box_z`` are not
+provided or set to 0, no periodic boundaries are applied