Includes various ways to calculate point spread functions (PSFs) with a focus on high NA fluorescence microscopy.
To run this toolbox, the user needs to install DipImage.
The MATLAB m files for generating the PSF are stored in the folder "Code". To run the code, change the current directory to the directory of the toolbox. This can be done using the code below:
cd <PSFToolbox-Matlab directory> If necessary, add the path of each subfolder in the Code folder into the MATLAB current directory.
Use the test.m file to run a test of the computation.
The main function for computing the PSF is called GenericPSFSim. It takes the following as input and output argument.
ImageParam: a structure containing the following fields:
- Sampling : voxel size represented in a vector format of three elements along x, y and, z
- Size : window size along x, y and, z
PSFParam: a structure containing the following minimum fields:
- polarization : polarization state of the field. It can take as an input integers 'linearX' or 'linearY' or 'radial' or 'azimuthal' or 'circular'. The default polarization stated is assumed to be 'circular' if not mentionned
- Aplanatic : value can be
$\beta = -1, 0$ , or 1
This corresponds to the aplanatic correction applied in the pupil plane. For the correction, the pupil plane is multiplied by$cos^{\beta/2}(\theta)$ ,$\beta$ being the input argument in the Aplanatic and$\theta$ is the opening or elevation angle of the beam. A value of$\beta = 0$ means there is no aplanatic correction applied onto the field. By default,$\beta = -1$
- NA : numerical aperture of the objective lens
- n : refractive index of the immersion medium
- lambdaEm : emission wavelength
Method: a string for selecting the method to use in the simulation. Choices are:
- 'RichardsWolff' (default): Using the method from the Richards & Wolff paper to generate the electric field as well as the image intensity [1, 2].
- 'RichardsWolffInt' : A fast method, if only the intensity is needed [1, 2].
- 'RWFast': Reworked Richards & Wolf bast on an chirp-Z transform along kz [3, 4]. This method is still under revision.
- 'SlicePropagation' : uses the angular spectrum method to propagate the field from the in-plane amplitude transfer function (ATF) using a fast Fourier transform (FFT). A jinc function is used to simulate a smooth aperture to reduce Fourier artefacts [4].
- 'CZT' : a chirped-Z transform-based (== Zoomed FFT) method to avoid the out-of-focus wrap-around problems of SlicePropagation [3, 4].
- 'VolumeShell' : this creates a part of the McCutchen Pupil directly in Fourier space using pre-computed Fourier-space interpolators which is pretty good in the central region of the PSF but (intentionally) degrading near the edge of the volume [4].
- 'SincR' : this creates the Fourier shell by a 3D FFT of a sinc(abs(R))
AddParams: a structure with the details on the refractive indices and thickeness of each layer in the system
- 'ns': refractive index (ri) of the sample. If the sample is a stratified medium, value attributed to this variable can be in a vector format listing the refractive index of each layer such that the first element is the ri of where the emitter is at i.e. the the ri of the layer furtherst away from the coverslip. The last element of the vector is the ri of the medium closest to the coverslip.
- 'ng': ri of the coverslip in real condition.
- 'ni': ri of the immersion medium in real condition.
- 'ni0': ri of the immersion medium in design condition.
- 'ts': thickness of the sample or position where the emitter is at, this should be with the same length as the vector of ns and within the same order as ns.
- 'tg': thickness of the coverslip in real condition.
- 'tg0': thickness of the coverslip in design condition.
- 'wd': working distance which is supposed to be the thickness of the immersion medium in design condition.
The Sampling, thicknesses and wavelengths must be in the same units.
If AddParams is empty i.e. AddParams == [] or not added to the function, the code will run as with the assumption that there is no interfaces with different refractive index mismatches yielding to aberrations into the system. In this case, the refractive index that is input in the dictionary PSFParam will be used troughout the calculation as the ri of the medium in which light is propagating.
Example of usage: AddParams = struct('ns', 1.33, 'ng', 1.518, 'ni', 1.516, 'ni0', 1.518, 'ng0', 1.518, 'ts', 2e3, 'tg', 1.7e5, 'tg0', 1.7e5, 'wd', 1.5e5);
AddPhase : additional 2D extra phase that can be added into the function.
h: resulting intensity PSF
amp4d : 4D complex amplitude distribution (if supported).
ImageParam = struct('Sampling', [40 40 50], 'Size', [128 128 32]);
PSFParam = struct('NA', 1.2, 'n', 1.33, 'lambdaEm', 520);
AddParams = struct('ns', 1.33, 'ng', 1.518, 'ni', 1.518, 'ng0', 1.518, 'ni0', 1.518, 'ts', 2e3, 'tg', 1.7e5, 'tg0', 1.7e5, 'wd', 1.5e5);
[h,amp4d] = GenericPSFSim(ImageParam, PSFParam, 'VolumeShell', AddParams);[1] Wolf E. Electromagnetic diffraction in optical systems-I. An integral representation of the image field. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences. 1959 Dec 15;253(1274):349-57.
[2] Richards B, Wolf E. Electromagnetic diffraction in optical systems, II. Structure of the image field in an aplanatic system. Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences. 1959 Dec 15;253(1274):358-79.
[3] Rabiner L, Schafer RW, Rader C. The chirp z-transform algorithm. IEEE transactions on audio and electroacoustics. 1969 Jun;17(2):86-92.
[4] Miora RH, Rohwer E, Kielhorn M, Sheppard CJ, Bosman G, Heintzmann R. Calculating Point Spread Functions: Methods, Pitfalls and Solutions. arXiv preprint arXiv:2301.13515v2. 2024 March 10.