MicroscPSF / MicroscPSF-Cpp

Fast and Accurate 3D PSF Computation for Fluorescence Microscopy
MIT License
0 stars 0 forks source link

Singular matrix when `NA = 0.8` and `ni = 1.0` #11

Open antonysigma opened 2 hours ago

antonysigma commented 2 hours ago

When a lower NA objective lens is specified, the adaptive linear solver implementation [1] reports the following warning:

// Matlab: Ci = J \ phase
auto Ci = solve(J.t(), phase.t());

warning: solve(): system is singular (rcond: 1.61799e-16); attempting approx solution

Microscope parameters:

struct microscope_params_t {
    Micron ti0 = 1200.0_um;     //!< Expected immersion medium thickness
    float ni0 = 1.0;           //!< Expected immersion medium refractive index
    float ni = 1.0;            //!< Measured immersion medium refractive index
    Micron tg0 = 170.0_um;     //!< Expected coverglass thickness
    Micron tg = 170.0_um;      //!< Measured coverglass thickness
    float ng0 = 1.5;           //!< Expected coverglass refractive index
    float ng = 1.5;            //!< Measured coverglass refractive index
    float ns = 1.33;           //!< Sample refractive index
    float NA = 0.8;            //!< Numerical aperture

    Micron pz = 2.0_um;  //!< Particle z distance from the sample-to-coverglass interface
};

auto voxel = scale_t<Micron>{.xy = 0.15, .z = 0.3};
auto volume = scale_t<uint32_t>{.xy = 128, .z = 63};
auto wavelength = 0.530_um;

[1] Sanderson and Curtin 2020, an adaptive solver for systems of linear equations. https://arma.sourceforge.net/armadillo_solver_2020.pdf

antonysigma commented 2 hours ago

Note to self: I notice a slight variation between the Matlab code and the Python code. Not sure if they are relevant to the singular matrix issue.

Currently, the scaling factor is computed as

const vec scaling_factor = (iota(1.0, precision.num_basis + 1) * 3 - 2) * params.NA *
                               (precision.min_wavelength / wavelength);

which follows the original literature and the code in MicroscPSF-py python code. The Matlab code is slightly different: the NA is normalized by the max NA parameter as params.NA / maxNA.

https://github.com/MicroscPSF/MicroscPSF-Matlab/blob/0ce640dbab71fdf4f0e7267091189e2a09e239d6/MicroscPSF.m#L145-L151