search

mehta-lab / waveorder

Wave optical models and inverse algorithms for label-agnostic imaging of density & orientation.
BSD 3-Clause "New" or "Revised" License
15 stars 4 forks source link

Odd-sized dimensions use incorrect frequency coordinates #90

Closed talonchandler closed 1 year ago

talonchandler commented 1 year ago

Spatial frequency coordinates are used by all of waveorder's frequency-space reconstructions (phase, uPTI, in-progress phase from polarization, fluorescence) to calculate the coordinates of OTFs and interpret the results of FTs.

For example, waveorder.util.gen_coordinate is used to calculate all transverse spatial frequency coordinates. If M is the number of pixels along an axis and ps is the pixel size, then the current implementation uses ifftshift((np.r_[:M]-M/2)/M/ps) to calculate the spatial frequencies along that dimension---see full function definition.

np.fft provides a utility function for exactly this purpose: fftfreq. The result of ifftshift((np.r_[:M]-M/2)/M/ps) and np.fft.fftfreq(M, ps) are identical for even M, but they diverge for odd M. The fftfreq function matches the convention of most DFT functions: if the signal has an even or odd length, the DC term should always appear as the first entry.

For example:

>>> import numpy as np
>>> M = 4
>>> ps = 1
>>> np.fft.ifftshift((np.r_[:M]-M/2)/M/ps)
array([0., 0.25, -0.5, -0.25])
>>> np.fft.fftfreq(M, ps)
array([0., 0.25, -0.5, -0.25]) # identical
>>> M = 5
>>> np.fft.ifftshift((np.r_[:M]-M/2)/M/ps)
array([-0.1,  0.1,  0.3, -0.5, -0.3]) # no DC term
>>> np.fft.fftfreq(M, ps)
array([ 0. ,  0.2,  0.4, -0.4, -0.2]) # different! now includes the DC term

This issue only affects reconstructions with odd dimensions, and the issue gets smaller as the dimensions get larger. I suspect a fix will only make a noticeable difference for small transverse FOVs, so this fix is a currently a low priority since most of our reconstructions are on large FOVs. .

I'm investigated the axial frequency coordinates further because we frequently use odd and small dimensions. I could not find any evidence of incorrect axial frequency coordinates for 3D-phase and 3D-fluorescence reconstructions because none of the direct computations happen in the frequency domain. I found potentially concerning lines in gen_dyadic_Greens_tensor which is used to generate the 3D vector WOTF for PTI, but I'm not extremely concerned because I think most of those reconstructions were run with a large number of z slices.

Context and comments welcome.

talonchandler commented 1 year ago

In alg-dev I'm using now using fftfreq to handle frequency-space coordinates in the generate_radial_frequencies function (used by all phase and fluorescence simulations and reconstructions), so this issue is solved.