Merge pull request #643 from mperrin/detector_offset

Allow shifting a detector, in alternative to shifting a source

poppy/poppy_core.py

@@ -1171,9 +1171,12 @@ def _propagate_mft(self, det):
         _log.debug(msg)
         self.history.append(msg)
 
+        if det.offset is not None:
+            _log.debug('      offset= '+str( det.offset))
         _log.debug('      MFT method = ' + mft.centering)
 
-        self.wavefront = mft.perform(self.wavefront, det_fov_lam_d, det_calc_size_pixels)
+        self.wavefront = mft.perform(self.wavefront, det_fov_lam_d, det_calc_size_pixels,
+                                     offset=None if det.offset is None else det.offset * det._offset_sign)  # sign flip intentional, see note in Detector class
         _log.debug("     Result wavefront: at={0} shape={1} ".format(
             self.location, str(self.shape)))
         self._last_transform_type = 'MFT'
@@ -3379,18 +3382,46 @@ class Detector(OpticalElement):
     oversample : int
         Oversampling factor beyond the detector pixel scale. The returned array will
         have sampling that much finer than the specified pixelscale.
+    offset : 2-tuple of floats
+        Offset (Y,X) in *pixels* for shifting the detector relative to the notional center of the output beam.
+        This has similar effect to shifting the source, but with opposite sign.
+        In other words, shifting a light source +1 arcsec in Y should have the same effect as
+        shifting the detector -1 arcsec in Y.
     """
 
     # Note, pixelscale argument is intentionally not included in the quantity_input decorator; that is
     # specially handled. See the _handle_pixelscale_units_flexibly method
     @utils.quantity_input(fov_pixels=u.pixel, fov_arcsec=u.arcsec)
     def __init__(self, pixelscale=1 * (u.arcsec / u.pixel), fov_pixels=None, fov_arcsec=None, oversample=1,
                  name="Detector",
+                 offset=None,
                  **kwargs):
         OpticalElement.__init__(self, name=name, planetype=PlaneType.detector, **kwargs)
         self.pixelscale = self._handle_pixelscale_units_flexibly(pixelscale, fov_pixels)
         self.oversample = oversample
 
+        if offset is not None:
+            if len(offset) != 2:
+                raise ValueError("If a detector offset is specified, it must be a tuple or list with 2 elements, "
+                                 "giving the (X, Y) offsets.")
+            # The offset is specified in pixels, so this can have units of pixels,
+            # or else if an integer or float, that's considered as implicitly a number of pixels
+            if isinstance(offset, u.Quantity):
+                try:
+                    offset = offset.to_value(u.pixel)
+                except u.UnitConversionError:
+                    raise(ValueError(f"A detector offset must be specified in units of detector pixels, not '{offset.unit}'"))
+            offset = np.asarray(offset)  # ensure it's an ndarray, not just a list or tuple
+        # A note on sign convention for detector offset: (This is regrettably confusing.)
+        #    The implementation in matrixDFT has the sense of "how much should the source be offset",
+        #    i.e. an offset of +5 pix moves the source by +5 pix.
+        #    However, physically we would like the opposite sign convention:  Moving the detector by +5 pix
+        #    should move the source by -5 pix.  This is implemented by a sign flip multplication by -1
+        #    which is applied in the _propagate_mft methods. That could just be a hard-coded -1,
+        #    but we choose to implement as a named variable to help make this logic clear later to readers of this code:
+        self.offset = offset
+        self._offset_sign = -1
+
         if fov_pixels is None and fov_arcsec is None:
             raise ValueError("Either fov_pixels or fov_arcsec must be specified!")
         elif fov_pixels is not None:

poppy/tests/test_core.py

@@ -710,6 +710,63 @@ def test_Detector_pixelscale_units():
             "Error message not as expected"
 
 
+def test_detector_offsets(plot=False, pixscale=0.01, fov_pixels=100):
+    """Test offsets of a detector.
+
+    It should be the case that:
+    (a) Offsettting the detector shifts the PSF
+    (b) And it does so with an opposite vector to shifting the source.
+    In other words, shifting a source by (+dX,+dY) should look the same as
+    shifting the detector by (-dX, -dY)
+
+    And you can specify the detector offsets in units of pixels or just as floats.
+
+    """
+    source_offset_r = .1
+    for with_units in [True, False]:
+        for offset_theta in [0, 45, 90, 180]:
+
+            # Compute offsets from radial to cartesian coords.
+            # recall astronomy convention is PA=0 is +Y, increasing CCW
+            source_offset_x = -source_offset_r * np.sin(np.deg2rad(offset_theta))  # arcsec
+            source_offset_y = source_offset_r * np.cos(np.deg2rad(offset_theta))
+            print(f"offset theta {offset_theta} is x = {source_offset_x}, y = {source_offset_y}")
+
+            # Create a PSF with a shifted source
+            offset_source_sys = poppy_core.OpticalSystem(npix=1024, oversample=1)
+            offset_source_sys.add_pupil(optics.ParityTestAperture())
+            offset_source_sys.add_detector(pixelscale=pixscale, fov_pixels=fov_pixels, oversample=1) #, offset=(pixscale/2, pixscale/2))
+            # This interface only has r, theta offsets available. Can't use _x, _y offsets here
+            offset_source_sys.source_offset_r = source_offset_r
+            offset_source_sys.source_offset_theta = offset_theta
+            offset_source_psf = offset_source_sys.calc_psf()
+
+            # Create a PSF with a shifted detector, the other way
+            offset_det_sys = poppy_core.OpticalSystem(npix=1024, oversample=1)
+            offset_det_sys.add_pupil(optics.ParityTestAperture())
+
+            det_offset = (-source_offset_y/pixscale, -source_offset_x/pixscale)   # Y, X in pixels
+            if with_units:
+                det_offset = np.asarray(det_offset) * u.pixel
+            offset_det_sys.add_detector(pixelscale=pixscale, fov_pixels=fov_pixels, oversample=1,
+                                        offset=det_offset)
+            offset_det_psf = offset_det_sys.calc_psf()
+
+            if plot:
+                fig, axes = plt.subplots(figsize=(16,9), ncols=3)
+                poppy.display_psf(offset_source_psf, ax=axes[0], crosshairs=True, colorbar=False,
+                                  title=f'Offset Source: {source_offset_x:.3f} arcsec, {source_offset_y:.3f}')
+                poppy.display_psf(offset_det_psf, ax=axes[1], crosshairs=True, colorbar=False,
+                                  title=f'Offset Det: {offset_det_sys.planes[-1].offset[1]:.2f},  {offset_det_sys.planes[-1].offset[0]:.2f} pix')
+                poppy.display_psf_difference(offset_source_psf, offset_det_psf, ax=axes[2],
+                                             title='difference', colorbar=False)
+
+            # Check the equality of the two results from the two above
+            assert np.allclose(offset_source_psf[0].data, offset_det_psf[0].data), "Offset source and offset detector the opposite way should be equivalent"
+
+
+
+
 # Tests for CompoundOpticalSystem