From 0078093dbd809128bef7fb62106cb13209bf0d5c Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Fri, 6 Dec 2024 16:36:45 +0100 Subject: [PATCH 1/4] Renormalize output array --- .../volume_rendering/off_axis_projection.py | 12 +++++++++--- 1 file changed, 9 insertions(+), 3 deletions(-) diff --git a/yt/visualization/volume_rendering/off_axis_projection.py b/yt/visualization/volume_rendering/off_axis_projection.py index 38c5677ef5..80ad9beefb 100644 --- a/yt/visualization/volume_rendering/off_axis_projection.py +++ b/yt/visualization/volume_rendering/off_axis_projection.py @@ -443,7 +443,7 @@ 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" ) @@ -451,12 +451,14 @@ def temp_weightfield(field, data): 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), From 6abd8fc72ff02d29e98c0dc6623bacbe06a25627 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Fri, 6 Dec 2024 17:09:54 +0100 Subject: [PATCH 2/4] Test w and w/out weighting --- .../tests/test_offaxisprojection.py | 54 +++++++++++-------- 1 file changed, 31 insertions(+), 23 deletions(-) diff --git a/yt/visualization/tests/test_offaxisprojection.py b/yt/visualization/tests/test_offaxisprojection.py index 46dccf860a..f72e3613ec 100644 --- a/yt/visualization/tests/test_offaxisprojection.py +++ b/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 bis + 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(): From 1ab751be0369a62ae212c92576e723adc7b30ca5 Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Fri, 6 Dec 2024 17:11:02 +0100 Subject: [PATCH 3/4] Discard the temporary weight field --- yt/visualization/volume_rendering/off_axis_projection.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/yt/visualization/volume_rendering/off_axis_projection.py b/yt/visualization/volume_rendering/off_axis_projection.py index 80ad9beefb..5de55f2fbb 100644 --- a/yt/visualization/volume_rendering/off_axis_projection.py +++ b/yt/visualization/volume_rendering/off_axis_projection.py @@ -498,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 = [] @@ -520,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: From 18ccf28c019fd19a1b1f977c92d41e3fa064b62d Mon Sep 17 00:00:00 2001 From: Corentin Cadiou Date: Fri, 6 Dec 2024 18:23:15 +0100 Subject: [PATCH 4/4] Fix typo Co-authored-by: Chris Havlin --- yt/visualization/tests/test_offaxisprojection.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/yt/visualization/tests/test_offaxisprojection.py b/yt/visualization/tests/test_offaxisprojection.py index f72e3613ec..39a640f629 100644 --- a/yt/visualization/tests/test_offaxisprojection.py +++ b/yt/visualization/tests/test_offaxisprojection.py @@ -239,7 +239,7 @@ def test_off_axis_octree(): v1, v2 = p1.frb["gas", "density"], p2.frb["gas", "density"] diff = (v1 - v2) / (v1 + v2) * 2 - # Make sure the difference has a small bis + # Make sure the difference has a small bias assert np.mean(diff).max() < 1e-3 # 0.1% # Compute 10-90% percentile