Bugfix: Fix normalization factors for gridder / imager

Makes more sense to keep beam / image vis at same normalization.

Then, post FFT we can correct for both sample-weightings and
FFT-normalization in a single operation.

Author: timstaley

src/fastimgproto/gridder/gridder.py

@@ -23,7 +23,10 @@ def convolve_to_grid(kernel_func,
                      raise_bounds=True,
                      progress_bar=None):
     """
-    Grid visibilities, calculating the exact kernel distribution for each.
+    Grid visibilities using convolutional gridding.
+
+    Returns the **un-normalized** weighted visibilities; the
+    weights-renormalization factor can be calculated by summing the sample grid.
 
     If ``exact == True`` then exact gridding is used, i.e. the kernel is
     recalculated for each visibility, with precise sub-pixel offset according to
@@ -135,13 +138,6 @@ def convolve_to_grid(kernel_func,
         if progress_bar is not None:
             progress_bar.update(1)
 
-    # Finally, renormalise the visibilities by the sampled weights total
-    # NB, this is computationally expensive, and can be side-stepped
-    # by simply tracking the sum of 'good_vis' weights, then normalizing
-    # source-fluxes later on. But we do it here for completeness.
-    sample_grid_total = sampling_grid.sum()
-    if sample_grid_total != 0. :
-        vis_grid /= sample_grid_total
     return vis_grid, sampling_grid
 
 

src/fastimgproto/imager.py

@@ -83,6 +83,16 @@ def image_visibilities(vis,
                                              )
     image = fft_to_image_plane(vis_grid)
     beam = fft_to_image_plane(sample_grid)
+
+    total_sample_weight = sample_grid.sum()
+    # To calculate the normalization factor:
+    # We correct for the FFT scale factor of 1/image_size**2
+    # (cf https://docs.scipy.org/doc/numpy/reference/routines.fft.html#implementation-details)
+    # And then divide by the total sample weight to renormalize the visibilities.
+    if total_sample_weight!=0:
+        renormalization_factor = (image_size_int*image_size_int)/total_sample_weight
+        beam *= renormalization_factor
+        image *= renormalization_factor
     if normalize:
         beam_max = np.max(np.real(beam))
         beam /= beam_max

tests/test_gridder/test_exact_gridding.py

@@ -36,7 +36,8 @@ def test_single_pixel_overlap_pillbox():
          [0., 0., 0., 0., 0., 0., 0., 0., ],
          [0., 0., 0., 0., 0., 0., 0., 0., ]]
     )
-    assert (expected_sample_grid * vis_amplitude == vis_grid).all()
+    assert (expected_sample_grid * vis_amplitude * vis_weights.sum()
+            == vis_grid).all()
     assert (expected_sample_grid * vis_weights.sum() == sampling_grid).all()
 
 
@@ -184,7 +185,7 @@ def test_multiple_complex_vis():
                                                vis_weights=vis_weights,
                                                )
     # simplification true since weights are all 1:
-    assert vis_grid.sum() == vis.sum() / vis_weights.sum()
+    assert vis_grid.sum() == vis.sum()
 
     # Since uv is precisely on a sampling point, we'll get a
     # 3x3 pillbox
@@ -201,7 +202,7 @@ def test_multiple_complex_vis():
          [0., v, v, v, 0., 0., 0., 0.],
          [0., 0., 0., 0., 0., 0., 0., 0.]]
     )
-    assert (expected_sample_grid / sampling_grid.sum() == vis_grid).all()
+    assert (expected_sample_grid  == vis_grid).all()
     assert (expected_sample_grid == sampling_grid).all()
 
 
@@ -227,7 +228,7 @@ def test_nearby_complex_vis():
                                                vis_weights=vis_weights,
                                                )
     # simplification true since weights are all 1:
-    assert vis_grid.sum() == vis.sum() / vis_weights.sum()
+    assert vis_grid.sum() == vis.sum()
 
     # Since uv is precisely on a sampling point, we'll get a
     # 3x3 pillbox
@@ -244,8 +245,8 @@ def test_nearby_complex_vis():
          [0., v, v, v, 0., 0., 0., 0.],
          [0., 0., 0., 0., 0., 0., 0., 0.]]
     )
-
-    assert (expected_sample_grid / vis_weights.sum() == vis_grid).all()
+    # simplification true since weights are all 1:
+    assert (expected_sample_grid == vis_grid).all()
     assert (expected_sample_grid == sampling_grid).all()
 
 

tests/test_gridder/test_oversampled_gridding.py

@@ -186,4 +186,4 @@ def test_oversampled_gridding():
                                )
 
     # simplification true since weights are all 1:
-    assert vis_grid.sum() == vis.sum() / vis_weights.sum()
+    assert vis_grid.sum() == vis.sum()

tests/test_gridder/test_sample_weighting.py

@@ -85,6 +85,6 @@ def test_natural_weighting():
          [0., 0., 0., 0., 0., 0., 0., 0., ]]
     )
     assert (
-    expected_sample_locations * vis_weights.sum() == sampling_grid).all()
+    expected_sample_locations == sampling_grid/sampling_grid.sum()).all()
     assert ((expected_sample_locations * natural_weighted_sum) ==
-                vis_grid).all()
+                vis_grid / sampling_grid.sum()).all()