bug in Projection Plots for SPH datasets

[nastasha-w]

Bug report

Bug summary

To my understanding, the pixel values in projection plots should represent the integral of some field over the line of sight through a pixel center. That means that pixels in different grids, but with the same center, should produce the same integrated values. However, this is not the case for the FixedResolutionBuffer projections of SPH datasets. I suspect this is due to the handling of small SPH particles (relative to the grid spacing) in the backend code: their contribution to the line-of-sight integrals seem to be systematically overestimated.

I've proposed a fix, but I'd like someone else to confirm whether this problem is real and seems correctly diagnosed before I open a pull request.

Bug explanation

Here, I will explain the issue as I understand it from the backend code. I will demonstrate the problem in the next section.

The backend code at issue here is pixelize_sph_kernel_projection in pixelization_routines.pyx, starting on line 1121. This code loops over the input particles, determines the extent of the region each particle might contribute to, and adds the (kernel-integrated path length) * (selected particle field) to each pixel's total. In #4783, I discuss some background on the reasoning for this method; the pixel values should be the same as the $\sum_{j} dl_j \cdot A_j$, along a Ray, where the sum is over SPH particles $j$, $dl_j$ is the path length from e.g., #4783, and $A$ is the particle field to be projected. For the column density or surface density $N$ along one of those pencil beams, we then get

N \approx \sum_{j=1}^{N'} A_j \cdot (m_j / (\rho_j \cdot h_j^2)) \cdot  \left[ \int dr \, U \left( \sqrt{(b_i/h_i)^2 + r^2}\right) \right],

and

dl_j = (m_j / (\rho_j \cdot h_j^2)) \cdot  \left[ \int dr \, U \left( \sqrt{(b_j/h_j)^2 + r^2} \right) \right].

For a particle $j$, $h_j$ is its smoothing length, $m_j$ is its mass, $\rho_j$ is its density, and $b_j$ is its impact parameter: its minimum distance to the line of sight.

The main loop in pixelize_sph_kernel_projection seems to implement thie. However, for small SPH particles, the smoothing lengths are tweaked in a way that seems to be inconsistent with the 'pencil-beam integral' approach to determining surface densities, and systematically overestimates their contributions to the surface density. This is because the calculation of $dl_j$ in the code is split up into two parts: a factor $A_j \cdot m_j / (\rho_j \cdot h_j^2)$ which is straightforwardly calculated from the dataset fields, and the integral $\int U \left( \sqrt{(b_j/h_j)^2 + r^2} dr \right)$, which only depends on $b_j/h_j$, and is therefore calculated for a small number of values, then interpolated for the rest with the SPHKernelInterpolationTable class.

This adjustment which I think is causing trouble is made in lines 1221-1223:

# we set the smoothing length squared with lower limit of the pixel
h_j2 = fmax(hsml[j]*hsml[j], dx*dy)
ih_j2 = 1.0/h_j2

where h_j2 is the smoothing length used in the kernel function integral, and ih_j2 is its inverse. hsml is the smoothing length recorded in the dataset, and dx and dy are the pixel dimensions.

However, the first factor in each particle's pixel contributions is calculated as

prefactor_j = pmass[j] / pdens[j] / hsml[j]**2 * quantity_to_smooth[j]

using the dataset smoothing length value hsml[j]. On the other hand, the integral is evaluated as

                            q_ij2 = (posx_diff + posy_diff) * ih_j2
                            if q_ij2 >= 1:
                                continue

                            # see equation 32 of the SPLASH paper
                            # now we just use the kernel projection
                            local_buff[xi + yi*xsize] +=  prefactor_j * itab.interpolate(q_ij2)

where posx_diff is (x-coordinate of particle center - x-coordinate of pixel center)^2, and posy_diff is the analogous value for y-coordinates. The integral evaluation (interpolation) happens with itab.interpolate(q_ij2). (This should indeed be a squared input value if I'm not mistaken; see #4783.)

For sufficiently large SPH particles (compared to the pixel grid spacing), this issue does not arise as fmax(hsml[j]*hsml[j], dx*dy) will be equal to the squared dataset smoothing length. Here, the different parts of the particle surface density calculation use consistent smoothing lengths. However, for smaller particles, the integral factor is evaluated at a smaller normalized impact parameter (b / hsml), and will therefore be larger than it should be.

For the inclusion of particles in a pixel calculation (i.e., the assignment non-zero $dl_j$ to particle $j$ for a given pixel center), the situation is complicated by the use of a mix of (h_j2, ih_j2) and hsml[j] in different parts of the calculation. To evaluate whether particle $j$ even lies in the domain of the projection, line 1232

if (px + hsml[j] < x_min) or (px - hsml[j] > x_max): continue

(and the analogous line 1236 for the y position) uses the true smoothing length hsml[j]. The other variables are the lower x-value bound in the grid xmin, the upper bound xmax, and px is the particle position x-coordinate (possibly with the box period added or subtracted to wrap around periodic edges). Starting on line 1238, the true smoothing length is also used to determine the square region of pixels particle $j$ might contribute to:

                    # here we find the pixels which this particle contributes to
                    x0 = <np.int64_t> ((px - hsml[j] - x_min)*idx)
                    x1 = <np.int64_t> ((px + hsml[j] - x_min)*idx)
                    x0 = iclip(x0-1, 0, xsize)
                    x1 = iclip(x1+1, 0, xsize)

                    y0 = <np.int64_t> ((py - hsml[j] - y_min)*idy)
                    y1 = <np.int64_t> ((py + hsml[j] - y_min)*idy)
                    y0 = iclip(y0-1, 0, ysize)
                    y1 = iclip(y1+1, 0, ysize)

The iclip function returns its first value if it's between the bounds given as its second and third arguments, but returns the closest bound if it is not. (I'm assuming second arg < thrid arg here.) idx and idy are the inverse pixel lengths in the x- and y-directions, respectively.

Two nested loops then iterate: for xi in range(x0, x1) and for yi in range(y0, y1). Even for a smoothing length of zero, this will return a range of at least 2 pixels (xi and yi values) in each coordinate direction. These nested loops start on line 1249:

                    # found pixels we deposit on, loop through those pixels
                    for xi in range(x0, x1):
                        # we use the centre of the pixel to calculate contribution
                        x = (xi + 0.5) * dx + x_min

                        posx_diff = px - x
                        posx_diff = posx_diff * posx_diff

                        if posx_diff > h_j2: continue

                        for yi in range(y0, y1):
                            y = (yi + 0.5) * dy + y_min

                            posy_diff = py - y
                            posy_diff = posy_diff * posy_diff
                            if posy_diff > h_j2: continue

                            q_ij2 = (posx_diff + posy_diff) * ih_j2
                            if q_ij2 >= 1:
                                continue

                            # see equation 32 of the SPLASH paper
                            # now we just use the kernel projection
                            local_buff[xi + yi*xsize] +=  prefactor_j * itab.interpolate(q_ij2)
                            mask[xi, yi] = 1

Within these loops, the distance between the position of particle $j$ and the pixel center is consistently compared to ih_j2 or h_j2: the adjusted smoothing lengths. This means that all SPH particles in the grid domain will contribute a non-zero density to some line of sight within the grid, even if they don't actually intersect any of the sightlines in this grid. Small SPH particles that do intersect the grid will have a larger $dl_j$ value (due to the integral part of the calculation) than their true smoothing lengths would give.

This means that the grid projections will systematically overestimate the contribution of small SPH particles (compared to the pixel grid spacing) to the line-of-sight integrals through pixel centers.

Note that this issue will also affect off-axis SPH projections, since the backend for those (off_axis_projection_SPH, starting on line 1915 of pixelization_routines.pyx) calls pixelize_sph_kernel_projection` after rotating the particle coordinates.

Proposed solution

If I am correct in diagnosing this issue, the fix would be simple: don't set a minimum effective smoothing length, and just set

h_j2 = hsml[j]*hsml[j]
ih_j2 = 1.0/h_j2

instead.

This does mean that some small SPH particles will be entirely 'missed' by the projection calculation. I think this ok for a grid of sightliness: this method will only conserve mass in the projection in the limit of a very high grid density, but it does sample the range of surface/column densities values without 'smoothing out' the extreme values.

An alternative approach to projecting SPH data does conserve mass explicitly, at the cost of 'smoothing out' the extremes of the
surface/column density range. Here, instead of projecting densities onto lines of sight through pixel centers, the mass of each SPH particle is divided up into different rectangular prisms (the pixels of the grid, extending along the projection axis). Here, weights can calculated by evaluating the kernel at pixel centers, then normalizing the sum of these weights to 1 explicitly. Here, setting a minimum smoothing length is a straightforward way to make sure each particle contributes to at least one pixel, and its mass is not lost. However, I do not believe this 'loss avoidance' strategy is appropriate for the 'sightline through the center' approach to projections used in YT.

Code for reproduction

Unfortunately, setting up the test case means this isn't quite a minimal amount of code, but it does straightforwardly suggest a test for this issue. Basically, in the first part, I just set up a grid of SPH particles with constant spacing along the x, y, and z coordinate directions, and a smoothing length equal to half the spacing, so the particles just 'touch' each other.

I then set up FixedResolutionBuffer grids of increasing pixel resolution, keeping the centers of a few pixels fixed. The centers coincide precisely with SPH particle centers; I set it up that way for a different test. Using these increasingly dense grids, I then retrieve the gas surface density at 9 pixel centers that are at the same coordinate locations in the different grids.

import numpy as np
from yt.loaders import load_particles

############ this part just sets up a dataset with SPH particles on a grid ###############
def constantmass(i: int, j: int, k: int) -> float:
    return 1.

## TODO: move to testing.py if this ends up in the test suite
def fake_sph_refineable_grid_ds(
        hsml_factor: float = 1.0,
        nperside: int = 3,
        periodic: bool = True, 
        e1hat: np.ndarray[float] = np.array([1, 0, 0]),
        e2hat: np.ndarray[float] = np.array([0, 1, 0]),
        e3hat: np.ndarray[float] = np.array([0, 0, 1]),
        offsets: np.ndarray[float] = 0.5 * np.ones((3,), dtype=np.float64),
        massgenerator = constantmass,
        unitrho: float = 1.,
        bbox: np.ndarray | None = None, 
        ):
    """Returns an in-memory SPH dataset useful for testing

    This dataset should have `nperside`**3 particles with the particles
    arranged uniformly on a 3D grid. 
    The basis vectors `e1hat`, `e2hat`, and `e3hat` define this grid.
    (If these are not normalized to 1 or not orthogonal, the spacing or
    overlap between SPH particles will be affected, but this is 
    allowed.) 
    Along these basis vectors, `offsets` give the zero point.
    All particles will have smoothing regions with a radius of 
    `hsml_factor` * 0.5,
    masses of 1, and densities of 1, and zero velocity.
    `unitrho` defines the density for a particle with mass 1 ('g'),
    and the standard (uniform) grid hsml_factor.
    """

    npart = nperside**3

    pos = np.empty((npart, 3), dtype=np.float64)
    mass = np.empty((npart,), dtype=np.float64)
    for i in range(0, nperside):
        for j in range(0, nperside):
            for k in range(0, nperside):
                _pos = (offsets[0] + i) * e1hat \
                        + (offsets[1] + j) * e2hat \
                        + (offsets[2] + k) * e3hat
                ind = nperside**2 * i + nperside * j + k
                pos[ind, :] = _pos
                mass[ind] = massgenerator(i, j, k)
    rho = unitrho * mass

    data = {
        "particle_position_x": (np.copy(pos[:, 0]), "cm"),
        "particle_position_y": (np.copy(pos[:, 1]), "cm"),
        "particle_position_z": (np.copy(pos[:, 2]), "cm"),
        "particle_mass": (mass, "g"),
        "particle_velocity_x": (np.zeros(npart), "cm/s"),
        "particle_velocity_y": (np.zeros(npart), "cm/s"),
        "particle_velocity_z": (np.zeros(npart), "cm/s"),
        "smoothing_length": (np.ones(npart) * 0.5 * hsml_factor, "cm"),
        "density": (rho, "g/cm**3"),
    }

    if bbox is None:
        eps = 1e-3
        margin = (1. + eps) * hsml_factor
        bbox = np.array([[np.min(_pos[:, 0]) - margin,
                          np.max(_pos[:, 0]) + margin],
                         [np.min(_pos[:, 1]) - margin,
                          np.max(_pos[:, 1]) + margin],
                         [np.min(_pos[:, 2]) - margin,
                          np.max(_pos[:, 2]) + margin], 
                         ])

    ds = load_particles(data=data, length_unit=1.0, 
                        bbox=bbox, periodicity=(periodic,) * 3)
    ds.kernel_name = 'cubic'
    return ds

def _get_dataset(reflevel: int = 1):
    '''
    constant density particle grid, 
    with increasing particle sampling
    '''
    lenfact = (1./3.)**(reflevel - 1)
    massfact = lenfact**3
    nperside = 3**reflevel

    e1hat = np.array([1, 0, 0]) * lenfact
    e2hat = np.array([0, 1, 0]) * lenfact
    e3hat = np.array([0, 0, 1]) * lenfact
    hsml_factor = lenfact
    bbox = np.array([[0., 3.]] * 3)
    offsets = np.ones(3, dtype=np.float64) * 0.5 * lenfact

    def refmass(i: int, j: int, k: int) -> float:
        return massfact
    unitrho = 1. / massfact # want density 1 for decreasing mass

    ds = fake_sph_refineable_grid_ds(hsml_factor=hsml_factor, 
                                     nperside=nperside,
                                     periodic=True, 
                                     e1hat=e1hat,
                                     e2hat=e2hat, 
                                     e3hat=e3hat,
                                     offsets=offsets,
                                     massgenerator=refmass,
                                     unitrho=unitrho,
                                     bbox=bbox, 
                                    )
    return ds 
################################################################

# actual test running
def test_gridproj3():
    # refine the pixel grid instead of the particle grid
    ds = _get_dataset(reflevel=2)
    proj = ds.proj(("gas", "density"), 2) 
    imgs = {}
    for rl in range(1, 6):
        npix = 1 + 2**(rl + 1)
        margin = 0.5 - 0.5**(rl + 1)
        frb = proj.to_frb(width=(3. - 2. * margin, 'cm'), 
                          resolution=(npix, npix),
                          height=(3. - 2. * margin, 'cm'),
                          center=np.array([1.5, 1.5, 1.5]), 
                          periodic=False)
        out = frb.get_image(('gas', 'density'))
        imgs[rl] = out
    return imgs

def run_main():
    imgs = test_gridproj3()
    print('Using a constant datset with pixel grids of increasing size, '
          'then selecting pixels with the same centers:')
    for rl in imgs:
        img = imgs[rl]
        pixspace = 2**(rl)
        print(f'Grid refinement level {rl}:') 
        print(img[::pixspace, ::pixspace])

if __name__ == '__main__':
    run_main() 

Actual outcome

>>> python sph_pixelization_tests.py
Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.
Intel MKL WARNING: Support of Intel(R) Streaming SIMD Extensions 4.2 (Intel(R) SSE4.2) enabled only processors has been deprecated. Intel oneAPI Math Kernel Library 2025.0 will require Intel(R) Advanced Vector Extensions (Intel(R) AVX) instructions.
yt : [INFO     ] 2024-01-25 01:24:28,451 Parameters: current_time              = 0.0
yt : [INFO     ] 2024-01-25 01:24:28,451 Parameters: domain_dimensions         = [1 1 1]
yt : [INFO     ] 2024-01-25 01:24:28,452 Parameters: domain_left_edge          = [0. 0. 0.]
yt : [INFO     ] 2024-01-25 01:24:28,452 Parameters: domain_right_edge         = [3. 3. 3.]
yt : [INFO     ] 2024-01-25 01:24:28,452 Parameters: cosmological_simulation   = 0
yt : [INFO     ] 2024-01-25 01:24:28,453 Allocating for 729 particles
yt : [INFO     ] 2024-01-25 01:24:28,765 Making a fixed resolution buffer of (('gas', 'density')) 5 by 5
yt : [INFO     ] 2024-01-25 01:24:28,774 Making a fixed resolution buffer of (('gas', 'density')) 9 by 9
yt : [INFO     ] 2024-01-25 01:24:28,775 Making a fixed resolution buffer of (('gas', 'density')) 17 by 17
yt : [INFO     ] 2024-01-25 01:24:28,777 Making a fixed resolution buffer of (('gas', 'density')) 33 by 33
yt : [INFO     ] 2024-01-25 01:24:28,778 Making a fixed resolution buffer of (('gas', 'density')) 65 by 65
Using a constant datset with pixel grids of increasing size, then selecting pixels with the same centers:
Grid refinement level 1:
[[26.17550883 26.17550883 26.17550883]
 [26.17550883 26.17550883 26.17550883]
 [26.17550883 26.17550883 26.17550883]] g/cm**2
Grid refinement level 2:
[[2.05048465 2.05048465 2.05048465]
 [2.05048465 2.05048465 2.05048465]
 [2.05048465 2.05048465 2.05048465]] g/cm**2
Grid refinement level 3:
[[0.00345061 0.00345061 0.00345061]
 [0.00345061 0.00345061 0.00345061]
 [0.00345061 0.00345061 0.00345061]] g/cm**2
Grid refinement level 4:
[[0.00345061 0.00345061 0.00345061]
 [0.00345061 0.00345061 0.00345061]
 [0.00345061 0.00345061 0.00345061]] g/cm**2
Grid refinement level 5:
[[0.00345061 0.00345061 0.00345061]
 [0.00345061 0.00345061 0.00345061]
 [0.00345061 0.00345061 0.00345061]] g/cm**2

Expected outcome

My expectation was that a pixel value represents the integral of a quantity (e.g., gas density) over a line of sight (pencil beam) through the pixel center. This means that the positions of the other pixels should not affect the value attributed to a pixel through the same center, e.g., if that center shows up in different pixel grids like in my test.

However, those values are only constant when the grid spacing is smaller than the SPH particle sizes. At lower pixel grid resolution, the surface densities come out too high.

Version Information

• Operating System: Mac OS Sonoma 14.2.1 (23C71)
• Python Version: 3.11.7
• yt version: 4.4.dev0
• numpy version: 1.26.3

I installed most dependencies from conda (standard channel), but installed yt from source, using my conda environment's pip.

[nastasha-w]

I think basically the same issue also comes up in slices through SPH datasets: pixelization_routines.pyx, line 1465 onwards: pixelize_sph_kernel_slice. The general approach here seems to be to evaluate the SPH field $A$ at the center $\vec{x}$ of each slice pixel:

A(\vec{x}) \approx \sum_{j=1}^{N} A_j \cdot (m_j / \rho_j) \cdot W\left(|\vec{x} - \vec{x_j}|, h_j\right),

where $A_j$ is the value of that field for particle $j$, $m_j$ is its mass, $\rho_j$ is its density, $\vec{x_j}$ is its position, and $h_j$ is its smoothing length. For small particles with smoothing lengths $\lesssim$ the pixel size, the smoothing lengths are artificially stretched, so that the normalized kernel and prefactor parts of the contribution are calculated inconsistently, and small SPH particles add to the total even when they do not intersect the pixel center.

Finally, I would suggest that we really need the z values for this as well. Now, this method calculates the kernel contributions for all particles as if the particle center lies exactly in the plane of the slice, but really, the kernel should probably be evaluated at the 3D-distance to the pixel center, to be consistent with how SPH interpolation is done in a simulation.

[nastasha-w]

The 3D version of this grid deposition pixelization_routines.pyx, line 1607 onwards: pixelize_sph_kernel_arbitrary_grid seems to handle the 3D part just find, but again sets minimum smoothing lengths (for the kernel evaluation, but not the pre-factor). This is inconsistent with an evaluation of the SPH field at a voxel center.

[nastasha-w]

I think I have the ProjectionPlot issues fixed on this branch: https://github.com/nastasha-w/yt/tree/sph_proj_backend. This means both the fixes to the on-axis projection plots for SPH datasets, and the addition of an off-axis projection option for SPH datasets. This new option just works through ProjectionPlot and OffAxisProjectionPlot, just like it does for grid data, with the same arguments, so I don't think we really need any new documentation there.

Since this took me so long, actually merging this back into main yt will likely be a bit painful.

[nastasha-w]

I did put the tests I used in a script (attached):
sph_pixelization_tests_2024-07-03.py.zip I haven't put those into the main test suite yet, but that should hopefully not be too difficult.

There is still the question of what to do with SPH particles that are less than a kernel size away from the projection region in depth, but are centered outside the region. YT regions, rays, etc. include any SPH particle that intersects region, but projections integrate the kernel over the line of sight assuming the line of sight goes all the way through the kernel. That means that, basically, the line of sight if effectively a bit longer than specified. For a projection around e.g., a single halo, where the edges are typically in some low-density region, that's probably not a problem, but in projections through full periodic volumes, that likely means double-counting particles less than a smoothing length from the boundary in the line-of-sight direction. I think this can be fixed without breaking the projections for smaller regions overlapping the periodic edges by shifting the line-of-sight coordinates (and the boundaries) to the projection region center and then ignoring the smoothing length 'margin' along the line of sight direction.

[jzuhone]

@nastasha-w can you submit a PR for this? I think that this will be easier to look at this way. I'm ready and willing to start looking at it next week.

[nastasha-w]

Thanks! Here's the PR: #4939

[nastasha-w]

Based on some tests with PR #4939 (this comment), I can now confirm that on test datasets from OWLS and Arepo, the resolution effects I showed in a test case in this issue play out in real data analysis, i.e., the bug is no longer just 'suspected'.

[jzuhone]

Could I get a weigh-in from @JBorrow on this issue and PR #4939 ? You are far more expert on these issues than I.

[JBorrow]

Will take a look. Based on my initial quick look over the issue, this paper: https://ui.adsabs.harvard.edu/abs/2021arXiv210605281B/abstract will likely be relevant. More soon.

[nastasha-w]

Thanks @JBorrow! I had seen that paper a while ago. The backend here uses a bit of a different approach from that paper: instead of dividing the mass/number of ions/other integral quantity in each SPH particle across rectangular prisms defined by the pixel grid, the existing SPH projection back-end in yt integrates the density/ion density/per-unit-volume quantity over a pencil-beam line of sight through the pixel centers. I removed some parts of that code that were (I think wrongly) trying to 'correct' for particles small enough to potentially not intersect any of these sightlines by, in some places, attributing a larger smoothing length to them than the simulation data listed. With this pencil-beam approach, the point is to sample specific sightlines, so 'undersampling' SPH kernels is allowed. This has the downside that this projection approach is not necessarily mass-conserving (though it should do so in the limit where the grid spacing is >> the smallest smoothing lengths), but the upside that at a low grid resolution, column/surface densities are not 'oversmoothed' (averaged over pixel scales) compared to pencil-beam observations. Certainly, there are cases (e.g., when modeling emission instead of absorption) where the pixel-averaging approach would be more suitable, but I've already spent a few weeks on this PR, so I'm going to leave that on the 'I may or may not get to it' pile. If you are interested in putting that algorithm into yt as a projection option, I'd certainly like to have them both available though!

[nastasha-w]

Honestly, a bigger worry for me is currently the periodic edge handling. The way this currently works is that, if a simulation is periodic in a particular coordinate direction, particles are effectively 'doubled', i.e., placed on both relevant sides of the grid. Perpendicular to the line of sight direction, this works just fine. However, this is also applied along the line-of-sight direction, where I think it might destructively come together with the line-of-sight integration approach. When making projections, integrating the density of a single SPH particle comes down to integrating the kernel across the line of sight, then multiplying that by the density and some normalization factors. Since the kernel only depends on (distance / smoothing length), we can factor out a power of the smoothing length and get an integral that only depends on (impact parameter / smoothing length). In the yt backend, the developer who actually wrote this projection code made clever use of that: they tabulated the dimensionless integral values as a function of (impact parameter / smoothing length), then used linear interpolation to speed up the integral calculation for each particle. The rub is that the integral only has this single-parameter dependence if we're integrating all the way through the kernel, which is what the code currently does for every SPH particle along the line of sight. In other words, if the particle is 'doubled' along the projection direction, it might actually be double-counted if we're projecting through the entire periodic volume along the line-of-sight direction. Simply ignoring periodicity in the line-of-sight direction would be a bad idea, since this would break functionality for thinner projections near periodic boundaries as well. However, along the line-of-sight direction, I think removing the smoothing-length margin to check whether a particle intersects the grid would be a reasonable fix. Perpendicular to the line of sight, we do not want to do this, since we'd miss particles that intersect some grid sightlines entirely. Integrating the kernel through the whole sphere even if the SPH particle 'sticks out' of the integration range along the line of sight means we're already overestimating some particles' contributions. Similarly excluding some SPH particles that are centered outside the volume but 'stick into' it might mean we're not making full use of the resolution at the ends of each sightline, but I don't think it should cause answers to be systematically off, like they are in that double-counting case. These same concerns apply to Ray objects in yt: their path lengths are calculated assuming the line of sight 'sticks out of the kernel' on both sides, but any SPH particle with a kernel that intersects the line segment defining the Ray is included in the particle selection.

[nastasha-w]

I've addressed the periodicity concern above in PR #4939 now (comment describing the effect).