plot PSF and MTF

hexSimProcessor.py

@@ -40,7 +40,8 @@ def _allocate_arrays(self):
         self._dx = self.pixelsize / self.magnification  # Sampling in image plane
         self._res = self.wavelength / (2 * self.NA)
         self._oversampling = self._res / self._dx
-        self._dk = self._oversampling / (self.N / 2)  # Sampling in frequency plane
+        self._dk = self._oversampling / (self.N / 2)  # Sampling in frequency plane. This is unitless, not sure what is means exactly
+        self._k_to_cycles_per_um = self.NA / self.wavelength
         self._kx = np.arange(-self._dk * self.N / 2, self._dk * self.N / 2, self._dk, dtype=np.single)
         [self._kx, self._ky] = np.meshgrid(self._kx, self._kx)
         self._dx2 = self._dx / 2
@@ -60,8 +61,8 @@ def _allocate_arrays(self):
         self._lastN = self.N
 
     def calibrate(self, img):
-        if self.N != self._lastN:
-            self._allocate_arrays()
+        # Always allocate, since not just N can change and thus res, kx, ... have to be adjusted
+        self._allocate_arrays()
 
         kr = sqrt(self._kx ** 2 + self._ky ** 2, dtype=np.single)
         kxbig = np.arange(-self._dk * self.N, self._dk * self.N, self._dk, dtype=np.single)
@@ -98,7 +99,9 @@ def calibrate(self, img):
             print(f'p  = {p[0]}, {p[1]}, {p[2]}')
             print(f'a  = {ampl[0]}, {ampl[1]}, {ampl[2]}')
 
-        ph = np.single(2 * pi * self.NA / self.wavelength)
+        ph = np.single(2 * pi * self.NA / self.wavelength)  # converts k to rad per um
+
+        self.mean_carrier_freq = np.mean(np.sqrt(ckx ** 2 + cky ** 2)) * ph  # [rad/um]
 
         xx = np.arange(-self._dx2 * self.N, self._dx2 * self.N, self._dx2, dtype=np.single)
         yy = xx
@@ -343,6 +346,7 @@ def _attm(self, kr, mask):
         return atff.reshape(kr.shape)
 
     def _tf(self, kr):
+        # TODO according to wikipedia it should be 2/pi. Then it would also be 1 at freq 0. Why is it 1/pi here?
         otf = (1 / pi * (arccos(kr / 2) - kr / 2 * sqrt(1 - kr ** 2 / 4)))
         return otf
 

imaging_method.py

@@ -102,9 +102,24 @@ def __init__(self):
         self.reconstruction_offset_y = QtWidgets.QLineEdit("0")
         layout.addRow("Reconstruction offset y", self.reconstruction_offset_y)
 
+        self.reconstruction_pixelsize = QtWidgets.QLineEdit("11")
+        layout.addRow("Pixelsize [um]", self.reconstruction_pixelsize)
+
+        self.reconstruction_magnification = QtWidgets.QLineEdit("1.5")
+        layout.addRow("Magnification", self.reconstruction_magnification)
+
+        self.reconstruction_na = QtWidgets.QLineEdit("0.02")
+        layout.addRow("NA", self.reconstruction_na)
+
+        self.reconstruction_wavelength = QtWidgets.QLineEdit("0.660")
+        layout.addRow("Wavelength [um]", self.reconstruction_wavelength)
+
         self.reconstruction_eta_txt = QtWidgets.QLineEdit("0.5")
         layout.addRow("Reconstruction eta", self.reconstruction_eta_txt)
 
+        self.reconstruction_mtf_data = QtWidgets.QTextEdit("")
+        layout.addRow("MTF data", self.reconstruction_mtf_data)
+
         self.use_filter_chb = QtWidgets.QCheckBox("Use frequency space filtering")
         layout.addRow("Filter", self.use_filter_chb)
 
@@ -124,13 +139,15 @@ def __init__(self):
         self.fft_plot = PlotWidget()
         self.carrier_plot = PlotWidget()
         self.bands_plot = PlotWidget()
+        self.psf_plot = PlotWidget()
         self.recon_plot = PlotWidget(None, num_clim = 2)
 
         self.debug_tabs = [
             ('Recorded Images', self.image_plot),
             ('FFT', self.fft_plot),
             ('Carrier', self.carrier_plot),
             ('Bands', self.bands_plot),
+            ('PSF shaping', self.psf_plot),
             ('Reconstructed Image', self.recon_plot),
         ]
 
@@ -220,11 +237,23 @@ def calibrate(self):
         offset_x = int(self.reconstruction_offset_x.text())
         offset_y = int(self.reconstruction_offset_y.text())
         eta = float(self.reconstruction_eta_txt.text())
+#TODO save reconstruction parameters
+#TODO also maybe set these not as part of reconstruction but recording
+        pixelsize = float(self.reconstruction_pixelsize.text())
+        magnification = float(self.reconstruction_magnification.text())
+        NA = float(self.reconstruction_na.text())
+        wavelength = float(self.reconstruction_wavelength.text())
+        mtf_data = self.reconstruction_mtf_data.toPlainText()
+
 
         assert self.frames.shape[0] == 7
         frames = self.frames[:, offset_y:offset_y + N, offset_x:offset_x + N]
 
         self.p.N = N
+        self.p.pixelsize = pixelsize
+        self.p.magnification = magnification
+        self.p.NA = NA
+        self.p.wavelength = wavelength
         self.p.eta = eta
         self.p.calibrate(frames)
 
@@ -248,6 +277,49 @@ def calibrate(self):
         self.bands_plot.connect_clim(ims)
         self.bands_plot.plot.draw()
 
+
+        # TODO could do all this in pattern calculation also, would make more sense
+        #Expected PSF
+        def gaussian(x, mu, sig):
+            return np.exp(-np.power(x - mu, 2.) / (2 * np.power(sig, 2.)))
+
+        def plot_psf(plt, k):
+            N = 9
+            x = np.linspace(-self.p._dx*N/2, self.p._dx*N/2, N * 10)  # [um]
+            sigma = self.p._res / 2.355 # FWHM to sigma
+
+            # TODO res is not FWHM, or is it?. And PSF is not really gaussian but a sinc
+            psf_orig = gaussian(x, 0, sigma)
+            plt.plot(x, psf_orig, label="emission psf")
+
+            mod = (1 + np.cos(k * x)) / 2
+            plt.plot(x, mod, '--', label="modulation")
+
+            plt.plot(x, psf_orig * mod, label="SIM psf")
+
+            plt.vlines(np.arange(-self.p._dx*N/2, self.p._dx*N/2, self.p._dx), 0, 0.1, 'r', label="pixel in sample plane")
+
+            plt.set_xlabel("x [um]")
+            plt.set_title("PSF Shaping")
+            plt.legend(loc="upper right")
+
+        def plot_otf(plt):
+            plt.plot(self.p._kx[0] * self.p._k_to_cycles_per_um, self.p._tf(abs(self.p._kx[0])) * 2, label="OTF");
+            nyq = 0.5 / self.p._dx
+            plt.vlines([-nyq, nyq], 0, 0.2, label="nyquist frequency")
+            if mtf_data != "":
+                mtf = np.fromstring("\n".join(mtf_data.split("\n")[1:]), sep='\t').reshape(-1, 2)
+                plt.plot(mtf[:, 0] / 1000, mtf[:, 1])
+            plt.set_ylim(0, None)
+            plt.set_xlabel("Frequency [cycles / um]")
+            plt.legend(loc="upper right")
+
+        self.psf_plot.plot.fig.clear()
+        axs = self.psf_plot.plot.fig.subplots(1, 2)
+        plot_psf(axs[0], self.p.mean_carrier_freq)
+        plot_otf(axs[1])
+        self.psf_plot.plot.draw()
+
     def reconstruct(self):
         if not self.sim_enabled_chb.isChecked():
             return None, None