Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[BUG] Fix off-axis rendering when center is not [0.5, 0.5, 0.5], fix periodicity #5059

Merged
merged 2 commits into from
Dec 2, 2024
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
8 changes: 7 additions & 1 deletion yt/visualization/plot_window.py
Original file line number Diff line number Diff line change
Expand Up @@ -26,6 +26,7 @@
validate_moment,
)
from yt.geometry.api import Geometry
from yt.geometry.oct_geometry_handler import OctreeIndex
from yt.units.unit_object import Unit # type: ignore
from yt.units.unit_registry import UnitParseError # type: ignore
from yt.units.yt_array import YTArray, YTQuantity
Expand Down Expand Up @@ -2494,7 +2495,12 @@ def __init__(
is_sph_field = finfo.is_sph_field
particle_datasets = (ParticleDataset, StreamParticlesDataset)

if isinstance(data_source.ds, particle_datasets) and is_sph_field:
dom_width = data_source.ds.domain_width
cubic_domain = dom_width.max() == dom_width.min()

if (isinstance(data_source.ds, particle_datasets) and is_sph_field) or (
isinstance(data_source.ds.index, OctreeIndex) and cubic_domain
):
center_use = parse_center_array(center, ds=data_source.ds, axis=None)
else:
center_use = center_rot
Expand Down
30 changes: 30 additions & 0 deletions yt/visualization/tests/test_offaxisprojection.py
Original file line number Diff line number Diff line change
Expand Up @@ -16,6 +16,7 @@
from yt.visualization.api import (
OffAxisProjectionPlot,
OffAxisSlicePlot,
ProjectionPlot,
)
from yt.visualization.image_writer import write_projection
from yt.visualization.volume_rendering.api import off_axis_projection
Expand Down Expand Up @@ -210,6 +211,35 @@ def test_field_cut_off_axis_octree():
assert_equal(np.nanmin(p4rho[p4rho > 0.0]) >= 0.5, True)


def test_off_axis_octree():
ds = fake_octree_ds()
p1 = ProjectionPlot(
ds,
"x",
("gas", "density"),
center=[0.6] * 3,
width=0.8,
weight_field=("gas", "density"),
)
p2 = OffAxisProjectionPlot(
ds,
[1, 0, 0],
("gas", "density"),
center=[0.6] * 3,
width=0.8,
weight_field=("gas", "density"),
)

# Note: due to our implementation, the off-axis projection will have a
# slightly blurred cell edges so we can't do an exact comparison
v1, v2 = p1.frb["gas", "density"], p2.frb["gas", "density"]
diff = (v1 - v2) / (v1 + v2) * 2
chrishavlin marked this conversation as resolved.
Show resolved Hide resolved

# Make sure the difference is zero-centered with a small standard deviation
assert np.mean(diff).max() < 1e-3 # 0.1%: very little bias
assert np.std(diff) < 0.05 # <2% error on average


def test_offaxis_moment():
ds = fake_random_ds(64)

Expand Down
22 changes: 12 additions & 10 deletions yt/visualization/volume_rendering/off_axis_projection.py
Original file line number Diff line number Diff line change
Expand Up @@ -444,21 +444,23 @@ def temp_weightfield(field, data):
# We need the width of the plot window in projected coordinates,
# i.e. we ignore the z-component
wmax = width[:2].max()

# Normalize the positions & dx so that they are in the range [-0.5, 0.5]
xyz = np.stack(
[
((data_source["index", k] - center[i]) / wmax).to("1").d
for i, k in enumerate("xyz")
],
axis=-1,
xyz = data_source.ds.arr(
np.zeros((len(data_source[vol.field]), 3)), "code_length"
)

for idim, periodic in enumerate(data_source.ds.periodicity):
axis = data_source.ds.coordinates.axis_order[idim]
# Recenter positions w.r.t. center of the plot window
xyz[..., idim] = data_source["index", axis] - center[idim]
if not periodic:
continue
# Wrap into [-0.5, +0.5]
xyz[..., idim] = (xyz[..., idim] + 0.5) % 1 - 0.5
# If we have periodic boundaries, we need to wrap the corresponding
# coordinates into [-w/2, +w/2]
w = data_source.ds.domain_width[idim]
xyz[..., idim] = (xyz[..., idim] + w / 2) % w - w / 2

# Rescale to [-0.5, +0.5]
xyz = (xyz / wmax).to("1").d

dx = (data_source["index", "dx"] / wmax).to("1").d

Expand Down
Loading