tlambert03 / PSFmodels

Python bindings for scalar and vectorial models of the 3D microscope point spread function.
GNU General Public License v3.0
44 stars 10 forks source link

Pixel oversampling when nxy is even produces asymmetric PSF #7

Open ptbrown1729 opened 1 year ago

ptbrown1729 commented 1 year ago

Description

I have a use case where it is helpful to work with even-sized PSF arrays. The documentation in PSFmodels says not to do this, but I thought it might be useful to look into what the problem is. I am expecting the output PSF's to always be symmetric about the center, which is at pixel (nxy-1)//2. I find that this is the case except when nxy is an even number.

It looks like part of the issue is in this code

  R = new double[npx_];
  int idx = 0;
  double xi, yi;
  for (int y = -xymax_; y <= xymax_; ++y)
  {
    for (int x = -xymax_; x <= xymax_; ++x)
    {
      xi = (double)x - xp_;
      yi = (double)y - yp_;
      R[idx] = sqrt(xi * xi + yi * yi);
      ++idx;
    }
  }
}

where xymax_ = ((nx_)*p.sf - 1) / 2;

If, for example nx_ = 4 and p.sf = 3 then the coordinates would get grouped (-5, -4, -3); (-2, -1, 0); (1, 2, 3); (4, 5, ?), and the asymmetry between the 0th and 2nd grouping is clear

It seems like this should work for the even case also if the loop went from

xymin_  = -((nx_ - 1) // 2) * p.sf - (p.sf//2) ;
xymax_ = ((nx_ + 2) // 2) * p.sf - (p.sf + 1)//2;

for the even case there should be p.sf more coordinate points which are greater than zero versus less than zero.

nx_=4 and p.sf = 3 would then get grouped (-4, -3, -2); (-1, 0, 1); (2, 3, 4); (5, 6, 7)

What I Did

import psfmodels
import matplotlib.pyplot as plt

dxy = 0.065
wvl = 0.488
na = 1.3
neven = 4
nodd = 3

# without oversampling
sf = 1
psf_even1 = psfmodels.vectorial_psf(0, neven, dxy, wvl=wvl, params={"NA": na, "sf": sf}, normalize=False)
psf_odd1 = psfmodels.vectorial_psf(0, nodd, dxy, wvl=wvl, params={"NA": na, "sf": sf}, normalize=False)

# with oversampling
sf = 3
psf_even3 = psfmodels.vectorial_psf(0, neven, dxy, wvl=wvl, params={"NA": na, "sf": sf}, normalize=False)
psf_odd3 = psfmodels.vectorial_psf(0, nodd, dxy, wvl=wvl, params={"NA": na, "sf": sf}, normalize=False)

figh = plt.figure()
ax = figh.add_subplot(2, 2, 1)
ax.imshow(psf_even1[0])
ax.set_title("PSF, even sf = 1")

ax = figh.add_subplot(2, 2, 2)
ax.imshow(psf_odd1[0])
ax.set_title("PSF, odd sf = 1")

ax = figh.add_subplot(2, 2, 3)
ax.imshow(psf_even3[0])
ax.set_title("PSF, even sf = 3")

ax = figh.add_subplot(2, 2, 4)
ax.imshow(psf_odd3[0])
ax.set_title("PSF, odd sf = 3")

Figure_1

tlambert03 commented 9 months ago

hey @ptbrown1729, I'm so sorry this one slipped through the cracks.

I suspect you've long moved on to something that works for you. You're absolutely correct that using even sized arrays wouldn't be a difficult change. that "warning" was largely inherited from some legacy code and I hadn't looked into recently.

was this issue essentially just a suggestion for modification? If so, thanks very much, I appreciate you looking into it and I can incorporate these changes

ptbrown1729 commented 9 months ago

@tlambert03 - no worries, and yes just a suggestion for modification/improvement. This is an issue that took me a long time to track down, and still seems like unexpected behavior.

Where I use these functions, I wrap them with code to throw an error if oversampling is requested with an even sized grid.

I guess there is a bigger issue with supporting even-sized PSF's at all, which is there is not a natural choice of center location. So the PSF starts to depend on what routine you use for convolving an image with it.