diff --git a/scopesim/effects/psfs/discrete.py b/scopesim/effects/psfs/discrete.py
index 5422c7ee..cc4eb48a 100644
--- a/scopesim/effects/psfs/discrete.py
+++ b/scopesim/effects/psfs/discrete.py
@@ -234,10 +234,10 @@ def apply_to(self, fov, **kwargs):
         if not fov.fields:
             return fov
 
-        if fov.image is None:
-            fov.image = fov.make_image_hdu().data
+        if fov.hdu is None or fov.hdu.data is None:
+            fov.hdu = fov.make_image_hdu()
 
-        old_shape = fov.image.shape
+        old_shape = fov.hdu.data.shape
 
         # Get kernels that cover this fov, and their respective masks.
         # Kernels and masks returned by .get_kernel as list of tuples.
@@ -252,7 +252,7 @@ def apply_to(self, fov, **kwargs):
                 kernel /= sum_kernel
 
             # image convolution
-            image = fov.image.astype(float)
+            image = fov.hdu.data.astype(float)
             kernel = kernel.astype(float)
             new_image = convolve(image, kernel, mode="same")
             if canvas is None:
@@ -267,7 +267,7 @@ def apply_to(self, fov, **kwargs):
 
         # reset WCS header info
         new_shape = canvas.shape
-        fov.image = canvas
+        fov.hdu.data = canvas
 
         # TODO: careful with which dimensions mean what
         if "CRPIX1" in fov.header:
diff --git a/scopesim/optics/fov.py b/scopesim/optics/fov.py
index 79224b54..464edfcc 100644
--- a/scopesim/optics/fov.py
+++ b/scopesim/optics/fov.py
@@ -87,9 +87,10 @@ def __init__(self, header, waverange, detector_header=None, cmds=None, **kwargs)
         self.fields = []
         self.spectra = {}
 
+        # These are apparently not supposed to be used?
         self.cube = None        # 3D array for IFU, long-lit, Slicer-MOS
-        self.image = None       # 2D array for Imagers
-        self.spectrum = None    # SourceSpectrum for Fibre-fed MOS
+        # self.image = None       # 2D array for Imagers
+        # self.spectrum = None    # SourceSpectrum for Fibre-fed MOS
 
         self._waverange = None
         self._wavelength = None
diff --git a/scopesim/tests/tests_effects/test_FVPSF.py b/scopesim/tests/tests_effects/test_FVPSF.py
index 113afa98..b55e2c9b 100644
--- a/scopesim/tests/tests_effects/test_FVPSF.py
+++ b/scopesim/tests/tests_effects/test_FVPSF.py
@@ -113,7 +113,8 @@ def test_convolution_with_central_psf_for_central_region(self, centre_fov,
         nax1, nax2 = centre_fov.header["NAXIS1"], centre_fov.header["NAXIS2"]
         centre_fov.hdu.data = np.zeros((nax2, nax1))
         centre_fov.hdu.data[::3, ::3] = 1
-        sum_orig = np.sum(centre_fov.hdu.data)
+        orig = centre_fov.hdu.data.copy()
+        sum_orig = np.sum(orig)
 
         fvpsf = FieldVaryingPSF(filename=str(mock_path / "test_FVPSF.fits"))
         fov_back = fvpsf.apply_to(centre_fov)
@@ -130,11 +131,17 @@ def test_convolution_with_fvpsfs_for_shifted_region(self, centre_fov,
         centre_fov.hdu.data = np.zeros((nax2, nax1))
         centre_fov.hdu.data[1::5, 1::5] = 1
         centre_fov.fields = [1]
-        sum_orig = np.sum(centre_fov.hdu.data)
+        orig = centre_fov.hdu.data.copy()
+        sum_orig = np.sum(orig)
 
         fvpsf = FieldVaryingPSF(filename=str(mock_path / "test_FVPSF.fits"))
         fov_back = fvpsf.apply_to(centre_fov)
 
+        # The FieldVaryingPSF was broken from early 2022 till late 2024
+        # because the pixels were not modified at all. This assert verifies
+        # that the data is not identical after applying the PSF, thereby
+        # preventing regressions.
+        assert not (orig == fov_back.hdu.data).all()
         assert np.sum(fov_back.hdu.data) == approx(sum_orig, rel=1E-2)
 
         if PLOTS:
@@ -151,7 +158,8 @@ def test_circular_fvpsf(self, basic_circular_fvpsf, mock_path):
         centre_fov.hdu.data[x, y] = 1
         # centre_fov.hdu.data[6:nax1-6:10, 6:nax1-6:10] = 1
         centre_fov.fields = [1]
-        sum_orig = np.sum(centre_fov.hdu.data)
+        orig = centre_fov.hdu.data.copy()
+        sum_orig = np.sum(orig)
 
         fvpsf = FieldVaryingPSF(
             filename=str(mock_path / "test_circular_fvpsf.fits"))
@@ -162,6 +170,7 @@ def test_circular_fvpsf(self, basic_circular_fvpsf, mock_path):
             plt.show()
 
         # print(np.sum(fov_back.hdu.data), sum_orig)
+        assert not (orig == fov_back.hdu.data).all()
         assert np.sum(fov_back.hdu.data) == approx(sum_orig, rel=1E-2)