Merge pull request #5077 from cphyc/fix/5071

[FIX] Correctly normalize off-axis projections for octree datasets

yt/visualization/tests/test_offaxisprojection.py

@@ -212,32 +212,40 @@ def test_field_cut_off_axis_octree():
 
 
 def test_off_axis_octree():
+    np.random.seed(12345)
     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"),
-    )
+    center = [0.4, 0.4, 0.4]
+
+    for weight in [("gas", "cell_mass"), None, ("index", "dx")]:
+        p1 = ProjectionPlot(
+            ds,
+            "x",
+            ("gas", "density"),
+            center=center,
+            width=0.8,
+            weight_field=weight,
+        )
+        p2 = OffAxisProjectionPlot(
+            ds,
+            [1, 0, 0],
+            ("gas", "density"),
+            center=center,
+            width=0.8,
+            weight_field=weight,
+        )
+
+        # 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
 
-    # 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
+        # Make sure the difference has a small bias
+        assert np.mean(diff).max() < 1e-3  # 0.1%
 
-    # 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
+        # Compute 10-90% percentile
+        q10, q90 = np.percentile(diff, q=(10, 90))
+        assert q10 > -0.02  # 2%: little up/down deviations
+        assert q90 < +0.02  # 2%: little up/down deviations
 
 
 def test_offaxis_moment():

yt/visualization/volume_rendering/off_axis_projection.py

@@ -443,20 +443,22 @@ def temp_weightfield(field, data):
         data_source.get_data(fields)
         # We need the width of the plot window in projected coordinates,
         # i.e. we ignore the z-component
-        wmax = width[:2].max()
+        wmax = width[:2].max().to("code_length")
         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]
+            xyz[..., idim] = (data_source["index", axis] - center[idim]).to(
+                "code_length"
+            )
             if not periodic:
                 continue
             # If we have periodic boundaries, we need to wrap the corresponding
             # coordinates into [-w/2, +w/2]
-            w = data_source.ds.domain_width[idim]
+            w = data_source.ds.domain_width[idim].to("code_length")
             xyz[..., idim] = (xyz[..., idim] + w / 2) % w - w / 2
 
         # Rescale to [-0.5, +0.5]
@@ -483,6 +485,10 @@ def temp_weightfield(field, data):
             Nx=resolution[0],
             Ny=resolution[1],
         )
+        # Note: since dx was divided by wmax, we need to rescale by it
+        projected_weighted_qty *= wmax.d / np.sqrt(3)
+        projected_weight *= wmax.d / np.sqrt(3)
+
         image = ImageArray(
             data_source.ds.arr(
                 np.stack([projected_weighted_qty, projected_weight], axis=-1),
@@ -492,6 +498,7 @@ def temp_weightfield(field, data):
             registry=data_source.ds.unit_registry,
             info={"imtype": "rendering"},
         )
+
     else:
         for grid, mask in data_source.blocks:
             data = []
@@ -514,12 +521,14 @@ def temp_weightfield(field, data):
             vol.sampler(pg, num_threads=num_threads)
 
         image = vol.finalize_image(camera, vol.sampler.aimage)
+
         image = ImageArray(
             image, funits, registry=data_source.ds.unit_registry, info=image.info
         )
 
-        if weight is not None:
-            data_source.ds.field_info.pop(("index", "temp_weightfield"))
+    # Remove the temporary weight field
+    if weight is not None:
+        data_source.ds.field_info.pop(("index", "temp_weightfield"))
 
     if method == "integrate":
         if weight is None: