Updating optics and scatterers for Tutorial Review. (#210)

* first_commit_review

* Fixing documentation for parameters and defining new parameters.

- Íscat is defined by setting the input/output polarization to circular and adding a phase_shift_correction term.
- Documentation for all variables

* Update optics.py

Author: SkariZ

deeptrack/holography.py

@@ -27,7 +27,7 @@ def get_propagation_matrix(shape, to_z, pixel_size, wavelength, dx=0, dy=0):
 
 
 class Rescale(Feature):
-    """Rescales an optical field by subtracting the real part of the field beofre multiplication.
+    """Rescales an optical field by subtracting the real part of the field before multiplication.
 
     Parameters
     ----------

deeptrack/optics.py

@@ -353,12 +353,13 @@ def _pupil(
 
         defocus = np.reshape(defocus, (-1, 1, 1))
         z_shift = defocus * np.expand_dims(z_shift, axis=0)
-
+        #print(self.pupil)
         if include_aberration:
             pupil = self.pupil
             if isinstance(pupil, Feature):
 
                 pupil_function = pupil(pupil_function)
+
             elif isinstance(pupil, np.ndarray):
                 pupil_function *= pupil
 
@@ -669,6 +670,12 @@ def get(self, illuminated_volume, limits, fields, **kwargs):
                 include_aberration=True,
                 **kwargs,
             )[0],
+            self._pupil(
+                volume.shape[:2],
+                defocus=[0],
+                include_aberration=True,
+                **kwargs,
+            )[0]
         ]
 
         pupil_step = np.fft.fftshift(pupils[0])
@@ -677,7 +684,7 @@ def get(self, illuminated_volume, limits, fields, **kwargs):
         light_in = self.illumination.resolve(light_in)
         light_in = np.fft.fft2(light_in)
 
-        K = 2 * np.pi / kwargs["wavelength"]
+        K = 2 * np.pi / kwargs["wavelength"]*kwargs["refractive_index_medium"]
 
         z = z_limits[1]
         for i, z in zip(index_iterator, z_iterator):
@@ -690,14 +697,17 @@ def get(self, illuminated_volume, limits, fields, **kwargs):
             light = np.fft.ifft2(light_in)
             light_out = light * np.exp(1j * ri_slice * voxel_size[-1] * K)
             light_in = np.fft.fft2(light_out)
-
-        shifted_pupil = np.fft.fftshift(pupils[-1])
+  
+        shifted_pupil = np.fft.fftshift(pupils[1])
         light_in_focus = light_in * shifted_pupil
-
+        #import matplotlib.pyplot as plt
+        #plt.imshow(light_in_focus.imag)
+        #plt.show()
         if len(fields) > 0:
             field = np.sum(fields, axis=0)
             light_in_focus += field[..., 0]
-
+        shifted_pupil = np.fft.fftshift(pupils[-1])
+        light_in_focus = light_in_focus * shifted_pupil
         # Mask to remove light outside the pupil.
         mask = np.abs(shifted_pupil) > 0
         light_in_focus = light_in_focus * mask
@@ -724,6 +734,136 @@ def get(self, illuminated_volume, limits, fields, **kwargs):
 Holography = Brightfield
 
 
+class ISCAT(Brightfield):
+    """Images coherently illuminated samples using ISCAT.
+
+    Images samples by creating a discretized volume, where each pixel
+    represents the effective refractive index of that pixel. Light is
+    propagated through the sample iteratively by first propagating the
+    light in the fourier space, followed by a refractive index correction
+    in the real space.
+
+    Parameters
+    ----------
+    illumination : Feature
+        Feature-set resolving the complex field entering the sample. Default
+        is a field with all values 1.
+    NA : float
+        The NA of the limiting aperature.
+    wavelength : float
+        The wavelength of the scattered light in meters.
+    magnification : float
+        The magnification of the optical system.
+    resolution : array_like[float (, float, float)]
+        The distance between pixels in the camera. A third value can be
+        included to define the resolution in the z-direction.
+    refractive_index_medium : float
+        The refractive index of the medium.
+    padding : array_like[int, int, int, int]
+        Pads the sample volume with zeros to avoid edge effects.
+    output_region : array_like[int, int, int, int]
+        The region of the image to output (x,y,width,height). Default
+        None returns entire image.
+    pupil : Feature
+        A feature-set resolving the pupil function at focus. The feature-set
+        receive an unaberrated pupil as input.
+    illumination_angle : float
+        The angle relative to the optical axis. Default is π radians in ISCAT.
+    amp_factor : float
+        The amplitude factor of the field. Default is 1. 
+        The relative amplitude off the illuminating field and the reference field.
+
+    """
+
+    def __init__(
+            self,
+            illumination_angle = np.pi,
+            amp_factor = 1, 
+            **kwargs
+            ):
+
+        super().__init__(
+            illumination_angle=illumination_angle,
+            amp_factor=amp_factor,
+            input_polarization="circular",
+            output_polarization="circular",
+            phase_shift_correction=True,
+            **kwargs
+            )
+        
+class Darkfield(Brightfield):
+    """Images coherently illuminated samples using Darkfield.
+
+    Images samples by creating a discretized volume, where each pixel
+    represents the effective refractive index of that pixel. Light is
+    propagated through the sample iteratively by first propagating the
+    light in the fourier space, followed by a refractive index correction
+    in the real space.
+
+    Parameters
+    ----------
+    illumination : Feature
+        Feature-set resolving the complex field entering the sample. Default
+        is a field with all values 1.
+    NA : float
+        The NA of the limiting aperature.
+    wavelength : float
+        The wavelength of the scattered light in meters.
+    magnification : float
+        The magnification of the optical system.
+    resolution : array_like[float (, float, float)]
+        The distance between pixels in the camera. A third value can be
+        included to define the resolution in the z-direction.
+    refractive_index_medium : float
+        The refractive index of the medium.
+    padding : array_like[int, int, int, int]
+        Pads the sample volume with zeros to avoid edge effects.
+    output_region : array_like[int, int, int, int]
+        The region of the image to output (x,y,width,height). Default
+        None returns entire image.
+    pupil : Feature
+        A feature-set resolving the pupil function at focus. The feature-set
+        receive an unaberrated pupil as input.
+    illumination_angle : float
+        The angle relative to the optical axis. Default is π/2 radians.
+
+    """
+
+    def __init__(
+            self, 
+            illumination_angle = np.pi/2, 
+            **kwargs
+            ):
+        super().__init__(
+            illumination_angle=illumination_angle,
+            **kwargs)
+
+    #Retrieve get as super
+    def get(self, illuminated_volume, limits, fields, **kwargs):
+        """Retrieve the darkfield image of the illuminated volume.
+
+        Parameters
+        ----------
+        illuminated_volume : array_like
+            The volume of the sample being illuminated.
+        limits : array_like
+            The spatial limits of the volume.
+        fields : array_like
+            The fields interacting with the sample.
+        **kwargs : dict
+            Additional parameters passed to the super class's get method.
+
+        Returns
+        -------
+        numpy.ndarray
+            The darkfield image obtained by calculating the squared absolute
+            difference from 1.
+        """
+
+        field = super().get(illuminated_volume, limits, fields, return_field=True, **kwargs)
+        return np.square(np.abs(field-1))
+
+
 class IlluminationGradient(Feature):
     """Adds a gradient in the illumination
 

deeptrack/scatterers.py

@@ -489,10 +489,12 @@ class MieScatterer(Scatterer):
     input_polarization: float or Quantity
         Defines the polarization angle of the input. For simulating circularly
         polarized light we recommend a coherent sum of two simulated fields. For
-        unpolarized light we recommend a incoherent sum of two simulated fields.
+        unpolarized light we recommend a incoherent sum of two simulated fields. 
+        If defined as "circular", the coefficients are set to 1/2.
     output_polarization: float or Quantity or None
         If None, the output light is not polarized. Otherwise defines the angle of the
         polarization filter after the sample. For off-axis, keep the same as input_polarization.
+        If defined as "circular", the coefficients are multiplied by 1. I.e. no change.
     L : int or str
         The number of terms used to evaluate the mie theory. If `"auto"`,
         it determines the number of terms automatically.
@@ -507,8 +509,16 @@ class MieScatterer(Scatterer):
         If True, the feature returns the fft of the field, rather than the
         field itself.
     coherence_length : float
-        The temporal coherence length of a partially coherent light given in meters. If None, the illumination is
-        assumed to be coherent.
+        The temporal coherence length of a partially coherent light given in meters. 
+        If None, the illumination is assumed to be coherent.
+    amp_factor : float
+        A factor that scales the amplification of the field. 
+        This is useful for scaling the field to the correct intensity. Default is 1.
+    phase_shift_correction : bool
+        If True, the feature applies a phase shift correction to the output field. 
+        This is necessary for ISCAT simulations. 
+        The correction depends on the k-vector and z according to the formula: 
+        arr*=np.exp(1j * k * z + 1j * np.pi / 2)
     """
 
     __gpu_compatible__ = True
@@ -540,6 +550,9 @@ def __init__(
         position_objective=(0, 0),
         return_fft=False,
         coherence_length=None,
+        illumination_angle=0,
+        amp_factor=1,
+        phase_shift_correction=False,
         **kwargs,
     ):
         if polarization_angle is not None:
@@ -569,6 +582,9 @@ def __init__(
             position_objective=position_objective,
             return_fft=return_fft,
             coherence_length=coherence_length,
+            illumination_angle=illumination_angle,
+            amp_factor=amp_factor,
+            phase_shift_correction=phase_shift_correction,
             **kwargs,
         )
 
@@ -617,7 +633,7 @@ def get_XY(self, shape, voxel_size):
     def get_detector_mask(self, X, Y, radius):
         return np.sqrt(X**2 + Y**2) < radius
 
-    def get_plane_in_polar_coords(self, shape, voxel_size, plane_position):
+    def get_plane_in_polar_coords(self, shape, voxel_size, plane_position, illumination_angle):
 
         X, Y = self.get_XY(shape, voxel_size)
         X = image.maybe_cupy(X)
@@ -633,9 +649,10 @@ def get_plane_in_polar_coords(self, shape, voxel_size, plane_position):
 
         # get the angles
         cos_theta = Z / R3
+        illumination_cos_theta=np.cos(np.arccos(cos_theta)+illumination_angle)
         phi = np.arctan2(Y, X)
 
-        return R3, cos_theta, phi
+        return R3, cos_theta, illumination_cos_theta, phi
 
     def get(
         self,
@@ -657,9 +674,11 @@ def get(
         return_fft,
         coherence_length,
         output_region,
+        illumination_angle,
+        amp_factor,
+        phase_shift_correction,
         **kwargs,
     ):
-
         # Get size of the output
         xSize, ySize = self.get_xy_size(output_region, padding)
         voxel_size = get_active_voxel_size()
@@ -683,9 +702,10 @@ def get(
         )
 
         # get field evaluation plane at offset_z
-        R3_field, cos_theta_field, phi_field = self.get_plane_in_polar_coords(
-            arr.shape, voxel_size, relative_position * ratio
+        R3_field, cos_theta_field, illumination_angle_field, phi_field = self.get_plane_in_polar_coords(
+            arr.shape, voxel_size, relative_position * ratio, illumination_angle
         )
+        
         cos_phi_field, sin_phi_field = np.cos(phi_field), np.sin(phi_field)
         # x and y position of a beam passing through field evaluation plane on the objective
         x_farfield = (
@@ -706,43 +726,55 @@ def get(
         cos_theta_field = cos_theta_field[pupil_mask]
         phi_field = phi_field[pupil_mask]
 
-        if isinstance(input_polarization, (float, int, Quantity)):
-
+        illumination_angle_field=illumination_angle_field[pupil_mask]
+        
+        if isinstance(input_polarization, (float, int, str, Quantity)):
             if isinstance(input_polarization, Quantity):
                 input_polarization = input_polarization.to("rad")
                 input_polarization = input_polarization.magnitude
 
-            S1_coef = np.sin(phi_field + input_polarization)
-            S2_coef = np.cos(phi_field + input_polarization)
+            if isinstance(input_polarization, (float, int)): 
+                S1_coef = np.sin(phi_field + input_polarization) 
+                S2_coef = np.cos(phi_field + input_polarization)
+
+            # If the input polarization is circular set the coefficients to 1/2.
+            elif isinstance(input_polarization, (str)):
+                if input_polarization == "circular":
+                    S1_coef = 1/2
+                    S2_coef = 1/2
 
         if isinstance(output_polarization, (float, int, Quantity)):
             if isinstance(input_polarization, Quantity):
                 output_polarization = output_polarization.to("rad")
                 output_polarization = output_polarization.magnitude
 
             S1_coef *= np.sin(phi_field + output_polarization)
-            S2_coef *= np.cos(phi_field + output_polarization)
+            S2_coef *= np.cos(phi_field + output_polarization) * illumination_angle_field
 
         # Wave vector
         k = 2 * np.pi / wavelength * refractive_index_medium
 
         # Harmonics
         A, B = coefficients(L)
-        PI, TAU = D.mie_harmonics(cos_theta_field, L)
+        PI, TAU = D.mie_harmonics(illumination_angle_field, L)
 
         # Normalization factor
         E = [(2 * i + 1) / (i * (i + 1)) for i in range(1, L + 1)]
 
         # Scattering terms
         S1 = sum([E[i] * A[i] * PI[i] + E[i] * B[i] * TAU[i] for i in range(0, L)])
         S2 = sum([E[i] * B[i] * PI[i] + E[i] * A[i] * TAU[i] for i in range(0, L)])
-
+        
         arr[pupil_mask] = (
             -1j
             / (k * R3_field)
             * np.exp(1j * k * R3_field)
             * (S2 * S2_coef + S1 * S1_coef)
-        )
+        ) / amp_factor
+        
+        # For phase shift correction (a multiplication of the field by exp(1j * k * z)).
+        if phase_shift_correction:
+            arr *= np.exp(1j * k * z + 1j * np.pi / 2)
 
         # For partially coherent illumination
         if coherence_length:
@@ -778,6 +810,7 @@ def get(
             ),
         )
         fourier_field = fourier_field * propagation_matrix * np.exp(-1j * k * offset_z)
+
         if return_fft:
             return fourier_field[..., np.newaxis]
         else: