Singular matrix when NA = 0.8 and ni = 1.0

[antonysigma]

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]

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

[hijizhou]

It should be fine either way, just a normalization in case the value is too large.

[antonysigma]

It should be fine either way, just a normalization in case the value is too large.

Thank you for the guidance. So I suppose the singular matrix warning with J = besselj<0>(scaling_factor * Rho); has nothing to do with scaling. I will study further.

[antonysigma]

Update: The issue is resolved by down-adjusting microscope_param_t::rho_samples from 2000 to 100. We need a lower value because of low NA. Closing the ticket.

The over-determined, low rank matrix J when rho_samples = 2000:

And the high-rank matrix J when rho_samples = 100:

---

C++ Patch to preview the matrix:

diff --git a/examples/generate-psf.cpp b/examples/generate-psf.cpp
index 77846c4..b9fff4a 100644
--- a/examples/generate-psf.cpp
+++ b/examples/generate-psf.cpp
@@ -11,19 +11,19 @@ int main() {
     using namespace ::units::literals;
 
     microscope_params_t params{};
-    params.NA = 1.4;
-    params.ti0 = 150.0_um;
-    params.ni = 1.5;
-    params.ni0 = 1.5;
-    params.pz = 2.0_um;
+    params.NA = 0.8;
+    params.ti0 = 900.0_um;
+    params.ni = 1.0;
+    params.ni0 = 1.0;
+    params.pz = 0.3_um;
 
     precision_li2017_t precision{};
-    precision.num_basis = 153;
-    precision.rho_samples = 1000;
+    precision.num_basis = 200;
+    precision.rho_samples = 2000;
 
     const auto psf =
-        makePSF(params, {0.1_um, 0.25_um}, {256, 128}, 0.610_um, precision);
-#ifdef  ARMA_USE_HDF5
+        makePSF(params, {0.1_um, 0.25_um}, {130, 63}, 0.530_um, precision);
+#if 0
     std::cout << "Saving volume to HDF5...\n";
     psf.save(hdf5_name("psf.h5", "psf", hdf5_opts::trans));
     std::cout << R"(Done.
diff --git a/microsc-psf/src/main.cpp b/microsc-psf/src/main.cpp
index 9387cbf..74674ed 100644
--- a/microsc-psf/src/main.cpp
+++ b/microsc-psf/src/main.cpp
@@ -137,7 +137,9 @@ makePSF(microscope_params_t params, scale_t<Micron> voxel, scale_t<uint32_t> vol
     {
         // Define the basis of Bessel functions.
         // Shape: number of basis function by number of rho samples.
-        auto&& J = besselj<0>(scaling_factor * Rho);
+        std::cout << "Max(scaling_factor * Rho) = " << max(scaling_factor) * max(Rho) << "\n";
+        const auto& J = besselj<0>(scaling_factor * Rho);
+        mat(J * 127/ max(abs(J.as_col())) + 127).save("matrix_J.pgm", pgm_binary);
 
         // Compute the approximation to the sampled pupil phase by finding the least squares
         // solution to the complex coefficients of the Fourier-Bessel expansion.