PhysicalOptics.jl

apply_rot_symmetry(arr, xs, xs_out, ys_out)

Takes a 1D array arr with rotational symmetry about the center point. arr should have values evaluates at the positions xs.

xs_out and ys_out are the positions for the output array. We interpolate arr and evaluate it at distances sqrt(xs_out[j]^2+ys_out[i]^2).

calc_NA(focal_length, diameter[, n])

Calculate the numerical aperture of a system with focal_length and diameter of the lens. Per default n=1 is the refractive index.

Examples

julia> calc_NA(100e-3, 200e-3)
1.0

julia> calc_NA(100e-3, 200e-3, 1.33)
1.33

calc_k(λ, n=1)

Calculate the value of the (angular) wave number with vacuum wavelength λ in medium with refractive index n. It holds: $ k = \kappa \cdot 2 \pi $

calc_κ(λ, n=1)

Calculate the value of the (non angular) wave number with vacuum wavelength λ in medium with refractive index n. It holds: $ \kappa = \frac{k}{2 \pi} $

calc_λ(k, n=1)

Calculate the vacuum wavelength λ from the angular wave number k in medium with refractive index n. It holds: $ \lambda = \frac{n \cdot 2 \pi}{k} $

center_extract(arr, new_size_array)

Extracts a center of an array. new_size_array must be list of sizes indicating the output size of each dimension. Centered means that a center frequency stays at the center position. Works for even and uneven. If length(new_size_array) < length(ndims(arr)) the remaining dimensions are untouched and copied.

Examples

julia> center_extract([1 2; 3 4], [1]) 
1×2 Array{Int64,2}:
 3  4
julia> center_extract([1 2; 3 4], [1, 1])
1×1 Array{Int64,2}:
 4
julia> center_extract([1 2 3; 3 4 5; 6 7 8], [2 2])
2×2 Array{Int64,2}:
 1  2
 3  4

center_pos(x)

Calculate the position of the center frequency. Size of the array is x

Examples

julia> center_pos(3)
2
julia> center_pos(4)
3

center_set!(arr_large, arr_small)

Puts the arr_small central into arr_large. The convention, where the center is, is the same as the definition as for FFT based centered. Function works both for even and uneven arrays.

Examples

julia> center_set!([1, 1, 1, 1, 1, 1], [5, 5, 5])
6-element Array{Int64,1}:
 1
 1
 5
 5
 5
 1

circ(arr, L, radius_aperture)

Apply circular aperture of radius radius_aperture to an array which has field width of L. circ! is also available.

julia> circ(ones((5, 5)), 5, 2.5)
5×5 Array{Float64,2}:
 0.0  0.0  1.0  0.0  0.0
 0.0  1.0  1.0  1.0  0.0
 1.0  1.0  1.0  1.0  1.0
 0.0  1.0  1.0  1.0  0.0
 0.0  0.0  1.0  0.0  0.0

circ(radius_aperture, x, y)

Respresenting an aperture. Being 0 if x^2 + y^2 > radius_aperture^2 and otherwise 1.

Examples

julia> circ(1, 0, 1)
1

julia> circ(1, 0.5, 0.5)
1

julia> circ(1, 1, 1)
0

circ(radius_aperture, r)

Respresenting an aperture. Being 0 if r^2 > radius_aperture^2 and otherwise 1.

Examples

julia> circ(1, 0)
1

julia> circ(1, 1)
1

julia> circ(1, 1.01)
0

conv(u, v[, dims]; shift=false, real_res=false)

Convolve u with v over dims dimensions.

Arguments

• u is an array in real space.
• v is the array to be convolved.
• Per default dims=1:ndims(v) means that we perform the convolution over all dimensions of v. If dims is an array with integers, we perform convolution only over these dimensions. Eg. dims=[1,3] would perform the convolution over the first and third dimension. Second dimension is not convolved.
• Per default shift=false therefore we assume that the center point of v is alreay ifftshifted to the first entry of the array.
• Per default real_res=false means that the output will be real. Otherwise we cut off the imaginary part.

Examples

1D with FFT over all dimensions. We choose v to be a delta peak. Therefore convolution should act as identity.

julia> u = [1 2 3 4 5]
1×5 Array{Int64,2}:
 1  2  3  4  5

julia> v = [0 0 1 0 0]
1×5 Array{Int64,2}:
 0  0  1  0  0

julia> conv(u, v)
1×5 Array{Complex{Float64},2}:
 4.0+0.0im  5.0+0.0im  1.0+0.0im  2.0+0.0im  3.0+0.0im

julia> conv(u, v, real_res=true) # since v is not ifftshifted with peak at the first entry, we see a wrong result.
1×5 Array{Float64,2}:
 4.0  5.0  1.0  2.0  3.0

julia> conv(u, v, shift=true, real_res=true)
1×5 Array{Float64,2}:
 1.0  2.0  3.0  4.0  5.0

julia> conv(u, ifftshift(v), real_res=true)
1×5 Array{Float64,2}:
 1.0  2.0  3.0  4.0  5.0

2D with FFT with different dims arguments.

julia> u = [1 2 3; 4 5 6]
2×3 Array{Int64,2}:
 1  2  3
 4  5  6

julia> v = [1 0 0; 1 0 0]
2×3 Array{Int64,2}:
 1  0  0
 1  0  0

julia> conv(u, v, [2])
2×3 Array{Complex{Float64},2}:
 1.0+0.0im  2.0+0.0im  3.0+0.0im
 4.0+0.0im  5.0+0.0im  6.0+0.0im

julia> conv(u, v, [2], real_res=true)
2×3 Array{Float64,2}:
 1.0  2.0  3.0
 4.0  5.0  6.0

julia> conv(u, v, [1, 2], real_res=true) # now we do a 2D convolution which is not identity anymore
2×3 Array{Float64,2}:
 5.0  7.0  9.0
 5.0  7.0  9.0

julia> conv(u, v, real_res=true) # same statement as above
2×3 Array{Float64,2}:
 5.0  7.0  9.0
 5.0  7.0  9.0

conv_otf(obj, otf [ , dims])

Performs a FFT-based convolution of an obj with otf. Wrapper for conv_v_ft(obj, otf, dims, real_res=true). Check the help of conv_v_ft for more details and examples.

conv_otf_r(obj, otf [, dims])

Performs a FFT-based convolution of an obj with an otf. Same arguments as conv_otf but with obj being real and otf=rfft(psf). All FFTs are computed with rfft and irfft.

conv_psf(obj, psf [, dims]; shift=false)

Convolve obj with psf over dims dimensions. Based on FFT convolution. This function calls conv, check the help of this method. Wrapper for conv(obj, psf, dims, shift=shift, real_res=true)

This convolution is based on the real space sliding scheme. The kernel is expected to have the center point in the way like FFT defines it

arr is the arr and kernel is the small kernel.

conv_v_ft(u, v_ft[, dims]; real_res=false)

Convolve u with v_ft over dims dimensions. Based on FFT convolution.

Arguments

• u is an array in real space.
• v_ft is the array to be convolved with in Fourier space. Therefore you have to check yourself that v was shifted correctly in real space.
• Per default dims=1:ndims(v_ft) means that we perform the convolution over all dimensions of v_ft. If dims is an array with integers, we perform convolution only over these dimensions. Eg. dims=[1,3] would perform the convolution over the first and third dimension. Second dimension is not convolved.
• Per default real_res=false means that the output will be real. Otherwise we cut off the imaginary part.

Examples

Convolution with delta peak is an identity operation

julia> u = [1 2 3 4 5]
1×5 Array{Int64,2}:
 1  2  3  4  5

julia> v = [1 0 0 0 0]
1×5 Array{Int64,2}:
 1  0  0  0  0

julia> v_ft = fft(v)
1×5 Array{Complex{Float64},2}:
 1.0+0.0im  1.0+0.0im  1.0+0.0im  1.0+0.0im  1.0+0.0im

julia> conv_v_ft(u, v_ft)
1×5 Array{Complex{Float64},2}:
 1.0+0.0im  2.0+0.0im  3.0+0.0im  4.0+0.0im  5.0+0.0im

julia> conv_v_ft(u, v_ft, real_res=true)
1×5 Array{Float64,2}:
 1.0  2.0  3.0  4.0  5.0

debye_focus(L, N, NA, z; λ=550e-9, n=1)

Evaluates the Debye integral to get the electrical field at position z with respect to the focus of the lens. z>0 means farther away from the lens than the focus and vice versa. Lens has a numerical aperture NA. λ is the wavelength and n the refractive index.

Returns a quadratic array with size (N, N). The side length is L, specifically fftpos(L, N) are the coordinates.

fft_view_tuple(x[, shift, dims])

fftpos(L, N)

Construct a range from -L/2 to L/2. However, we ensure that everything is centered around the center in a way that a FFT interpretes it correctly. For odd sequences it is indeed in the real center. For even sequences the center is at N/2 + 1.

Examples

julia> collect(fftpos(1, 4))
4-element Array{Float64,1}:
 -0.5
 -0.25
  0.0
  0.25

julia> collect(fftpos(1, 5))
5-element Array{Float64,1}:
 -0.5
 -0.25
  0.0
  0.25
  0.5

four_f_propagate(arr, L, f1, f2, NA; fft_plan=fft_plan, d=nothing)

Propagate the electrical field arr (field size L) from front focal plane to the back focal plane of a 4f system. The focal length of the first lens is f1 and the second lens f. Magnification of the system is then f2/f1. d is the distance from lens to the point of the arr in front of the lens. If d == nothing then we propagate from object plane to image plane.

fresnel_kernel(fx, fy, z, λ, n=1.0)

Calculates the Fresnel propagation kernel in Fourier space. fx, f_y are the spatial frequency. z is the propagation distance. λ is the wavelength and n the refractive index.

fresnel_number(a, L, λ=λ0)

Calculate the Fresnel number where a is the characteristic size, L the distance from screen to aperture and λ the wavelength.

get_index_for_pos(pos_arr, pos)

For a arr containing range like data, return the index which is closest to the pos value

Example

julia> pos = collect(fftpos(1, 5))
5-element Vector{Float64}:
 -0.5
 -0.25
  0.0
  0.25
  0.5

julia> get_index_for_pos(pos, 0.124)
3

julia> get_index_for_pos(pos, 0.1251)
4

julia> pos = collect(fftpos(1, 6))
6-element Vector{Float64}:
 -0.5
 -0.3333333333333333
 -0.16666666666666666
  2.2204460492503052e-17
  0.16666666666666669
  0.33333333333333337

julia> get_index_for_pos(pos, 0)
4

julia> get_index_for_pos(pos, 0.1)
5

get_indices_around_center(i_in, i_out)

A function which provides two output indices i1 and i2 where i2 - i1 = i_out The indices are chosen in a way that the set i1:i2 cuts the interval 1:i_in in a way that the center frequency stays at the center position. Works for both odd and even indices

jinc(x)

Computes the jinc function which is $\text{jinc} = \frac{J_1(x)}{x}$ where $J_1$ being the first Bessel function.

In future, probably switch to: https://github.com/JuliaMath/SpecialFunctions.jl/issues/264

Examples

julia> jinc(0.0)
0.5

julia> jinc(0.0im)
0.5 + 0.0im

julia> jinc(0.0001im)
0.5000000006250005 - 3.0616170016954073e-17im

julia> jinc(1.0f0)
0.44005057f0

julia> jinc.([0.0 0.5 1.0])
1×3 Array{Float64,2}:
 0.5  0.484537  0.440051

jinc_psf(psf_size, L, radius[, f]; λ=550e-9, shift=false)

Generate the normalized, incoherent 2D jinc PSF of a circular aperture. psf_size is output array shape. L is the width of the array expressed in meter. radius is the aperture radius in meter. Keyword arguments λ and f=100e-3 represent wavelength and focal length of the lens respectively.

Reference: Mertz, J. (2019). Introduction to Optical Microscopy (2nd ed.).

lens(f, L, size=(512, 512); λ=550e-9, n=1)

Return complex lens transmission value for a field of size and a length of L The lens has a focal length f and the surrounding medium has refractive index n. Wavelength is λ.

lens(f, x, y; λ=550e-9, n=1)

Return complex transmission value at position x and y for a lens. The lens has a focal length f and the surrounding medium has refractive index n. Wavelength is λ.

lens_propagate(arr, L, f; λ=550e-9, n=1, d=nothing)

Propagate an electrical field of width L at a distance d in front of a lens with focal length f in medium with refractive index n to the back focal plane of the lens. Per default d=nothing meaning that we set d=f. Assume same sampling and size in both dimensions. As return parameter we get the resulting field and the new field size. Essentially the lens (assuminig infinite large of that lens) performs a scaled Fourier transform. Therefore we get a new field size.

Reference

Based on:

• "Computational Fourier Optics. A MATLAB Tutorial", D. Voelz, (2011).
• Goodman, Joseph W. Introduction to Fourier optics

micro_lens_array(f, L, pitch, size=(512, 512); centered=true, λ=550e-9, n=1, aperture=true,
                 dtype=ComplexF32)

Return a micro lens array (MLA) with pitch size pitch in the field size of L. All lenslets have the same focal length f. As default centered=true means that the there is a micro lens centered around the center of the array. aperture=true indicates that there is a circular aperture cutting the field. Otherwise the aperture are square shaped of the same size as the pitch.

mla_propagate(E, L, fmla, p_mla, z1, z2)

Propagate an electrical field E placed a distance z1 in front of a microlns array, to a distance z2 behind the microlens array. L is the size of the field and fmla is the focal length of the lenslets. p_mla is microlenses pitch.

normabs2(arr)

Takes abs2 and normalizes maximum to 1

Examples

julia> normabs2([10.0, 12.1])
2-element Array{Float64,1}:
 0.6830134553650707
 1.0

plan_conv_r(psf [, dims])

Pre-plan an optimized convolution for array shaped like psf (based on pre-plan FFT) along the given dimenions dims. dims = 1:ndims(psf) per default. The 0 frequency of psf must be located at the first entry. We return first the otf (obtained by rfft(psf)). The second return is the convolution function pconv. pconv itself has two arguments. pconv(obj, otf) where obj is the object and otf the otf. This function achieves faster convolution than conv_psf(obj, psf).

Examples

julia> u = [1 2 3 4 5]
1×5 Array{Int64,2}:
 1  2  3  4  5

julia> v = [1 0 0 0 0]
1×5 Array{Int64,2}:
 1  0  0  0  0

julia> otf, pconv = plan_conv_r(v)
(Complex{Float64}[1.0 + 0.0im 1.0 + 0.0im … 1.0 + 0.0im 1.0 + 0.0im], PhysicalOptics.var"#conv#49"{FFTW.rFFTWPlan{Float64,-1,false,2,UnitRange{Int64}},AbstractFFTs.ScaledPlan{Complex{Float64},FFTW.rFFTWPlan{Complex{Float64},1,false,2,UnitRange{Int64}},Float64}}(FFTW real-to-complex plan for 1×5 array of Float64
(rdft2-rank>=2/1
  (rdft2-r2hc-rank0-x5)
  (dft-direct-5 "n1fv_5_avx2_128")), 0.2 * FFTW complex-to-real plan for 1×5 array of Complex{Float64}
(rdft2-rank>=2/1
  (rdft2-hc2r-rank0
    (rdft-rank0-iter-ci/1-x5))
  (dft-direct-5 "n1bv_5_avx2_128"))))

julia> pconv(u, otf)
1×5 Array{Float64,2}:
 1.0  2.0  3.0  4.0  5.0

point_source_propagate(L, size, point; λ=550e-9, n=1, dtype=ComplexF64)

Propagate a point source in a field of width L with array size size over a distance z in a medium with refractive index n. Point source is located at position point(x, y, z) and then propagated to a location z=0. z<0 in point means forward propagation.

This is based on the analytical solution in real space and not on Fourier space propagation. The latter one suffers from artifacts while microscopic large L. This function should be always preferred for point sources.

propagate!(arr, L, z; kernel=rs_kernel, λ=550e-9, n=1)

Propagate an electric field arr with a array size (in physical dimensions like meter etc) of L over the distance z. Per default kernel=rs_kernel (Rayleigh-Sommerfeld) is the propagation kernel. λ is the wavelength and n the refractive index.

quadratic(arr, L, diameter, [, diameter=(0,0)])

Apply a quadratic aperture of diameter to an array which has a field width of L. The aperture is shifted by Δx, Δy with respect to the center

There is also an in-place version quadratic!

julia> quadratic(ones((4, 4)), 10, 5)
4×4 Array{Float64,2}:
 0.0  0.0  0.0  0.0
 0.0  1.0  1.0  1.0
 0.0  1.0  1.0  1.0
 0.0  1.0  1.0  1.0

julia> quadratic(ones((5, 5)), 10, 5)
5×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  1.0  1.0  1.0  0.0
 0.0  1.0  1.0  1.0  0.0
 0.0  1.0  1.0  1.0  0.0
 0.0  0.0  0.0  0.0  0.0

julia> quadratic(ones((5, 5)), 5, 2.5, (0.0, 1.25)) 
5×5 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  1.0  1.0
 0.0  0.0  1.0  1.0  1.0
 0.0  0.0  1.0  1.0  1.0
 0.0  0.0  0.0  0.0  0.0

rayleigh_criterion(focal_length, diameter; λ=500e-9)

Calculates the resolution of a microscope according to Rayleigh. focal_length is the focal length of the lens. λ is the wavelength and diameter is the diameter of the aperture of the lens. See Wikipedia.

Examples

julia> rayleigh_criterion(100, 200)
3.05e-7

julia> rayleigh_criterion(100, 200, λ=100e-9)
6.099999999999999e-8

rayleigh_criterion(NA; λ=500e-9)

Calculates the resolution of a microscope according to Rayleigh. λ is the wavelength and NA the numerical aperture of the system. See Wikipedia.

Examples

julia> rayleigh_criterion(1.0)
3.05e-7

julia> rayleigh_criterion(1.0, λ=100e-9)
6.099999999999999e-8

rr(s; L=nothing)

Generate a array with values being the distance to the center pixel. s specifies the output size of the 2D array. L specifies the length. Can be a single number indicating the same length for both dimensions or a tuple specifying a length for each dimension separately.

Examples

julia> PhysicalOptics.rr((4,4))
4×4 Matrix{Float64}:
 2.82843  2.23607  2.0  2.23607
 2.23607  1.41421  1.0  1.41421
 2.0      1.0      0.0  1.0
 2.23607  1.41421  1.0  1.41421

julia> PhysicalOptics.rr((4,4), L=1)
4×4 Matrix{Float64}:
 0.707107  0.559017  0.5   0.559017
 0.559017  0.353553  0.25  0.353553
 0.5       0.25      0.0   0.25
 0.559017  0.353553  0.25  0.353553

julia> PhysicalOptics.rr((4,4), L=(5, 1))
4×4 Matrix{Float64}:
 2.54951  1.34629  0.5   1.34629
 2.51247  1.27475  0.25  1.27475
 2.5      1.25     0.0   1.25
 2.51247  1.27475  0.25  1.27475

julia> PhysicalOptics.rr((5,5), L=(1, 1))
5×5 Matrix{Float64}:
 0.707107  0.559017  0.5   0.559017  0.707107
 0.559017  0.353553  0.25  0.353553  0.559017
 0.5       0.25      0.0   0.25      0.5
 0.559017  0.353553  0.25  0.353553  0.559017
 0.707107  0.559017  0.5   0.559017  0.707107

rs_kernel(fx, fy, z, λ, n=1.0)

Calculates the Rayleigh-Sommerfeld propagation kernel in Fourier space. fx, f_y are the spatial frequencies. z is the propagation distance. λ is the wavelength and n the refractive index.

simple_psf(psf_size, radius; shift=false)

Generation of an approximate 2D PSF. psf_size is the output size of the PSF. The PSF will be centered around the point [1, 1], radius indicates the pupil diameter in pixel from which the PSF is generated.

Examples

julia> simple_psf([5, 5], 2)
5×5 Array{Float64,2}:
 0.36       0.104721    0.0152786    0.0152786    0.104721
 0.104721   0.0304627   0.00444444   0.00444444   0.0304627
 0.0152786  0.00444444  0.000648436  0.000648436  0.00444444
 0.0152786  0.00444444  0.000648436  0.000648436  0.00444444
 0.104721   0.0304627   0.00444444   0.00444444   0.0304627