Problem getting DECODE to work with an sCMOS camera (no EM gain)

[cleterrier]

Hi, I'm trying to train DECODE for acquired sequences from a microscope with sCMOS cameras rather than EMCCD.
I took care of changing the camera parameters in SMAP for the Z calibration and processing a sample of frames to obtain a Training parameter yaml file, which is the following (not that as I defined that there was no EM gain on the sCMOS camera, the em_gain value is 0.0):

Camera:
  baseline: 100.0
  e_per_adu: 0.48
  em_gain: 0.0
  px_size: [97.0, 97.0]
  read_sigma: 74.4
  spur_noise: 0.002
InOut: {calibration_file: 'D:/data SMLM/20210303 Calib 3D Karo/210209_2C_2Cameras_100nm/zStackCamera1_3dcal.mat',
  experiment_out: 'D:/data SMLM/20210303 Calib 3D Karo/210209_2C_2Cameras_100nm'}
Simulation:
  bg_uniform: [24.927, 38.2988]
  emitter_av: 25.0
  emitter_extent:
  - [-0.5, 39.5]
  - [-0.5, 39.5]
  - [-62.0, 507.0]
  intensity_mu_sig: [2284.5032, 456.90063]
  lifetime_avg: 1.0438

When using this yaml file as input in the training Jupyter notebook, first I get an error at the step that makes:

tar_em, sim_frames, bg_frames = simulator.sample()
sim_frames = sim_frames.cpu()

frame_path = 'D:/data SMLM/20210303 Calib 3D Karo/C1_MT-AF647_Clath-CF568_dual_astig_25ms_47um_27k_frames/ROI1_5K_crop40.tif'  # change if you load your own data
data_frames = decode.utils.frames_io.load_tif(frame_path).cpu()

print(f'Data shapes, simulation: {sim_frames.shape}, real data: {data_frames.shape}')
print(f'Average value, simulation: {sim_frames.mean().round()}, real data: {data_frames.mean().round()}')

The error being:

---------------------------------------------------------------------------
ZeroDivisionError                         Traceback (most recent call last)
<ipython-input-10-4d2b30743787> in <module>()
----> 1 tar_em, sim_frames, bg_frames = simulator.sample()
      2 sim_frames = sim_frames.cpu()
      3 
      4 frame_path = 'D:/data SMLM/20210303 Calib 3D Karo/C1_MT-AF647_Clath-CF568_dual_astig_25ms_47um_27k_frames/ROI1_5K_crop40.tif'  # change if you load your own data
      5 data_frames = decode.utils.frames_io.load_tif(frame_path).cpu()

C:\Users\chris\.conda\envs\decode_env\lib\site-packages\decode\simulation\simulator.py in sample(self)
     54 
     55         emitter = self.em_sampler()
---> 56         frames, bg = self.forward(emitter)
     57         return emitter, frames, bg
     58 

C:\Users\chris\.conda\envs\decode_env\lib\site-packages\decode\simulation\simulator.py in forward(self, em, ix_low, ix_high)
     93 
     94         if self.noise is not None:
---> 95             frames = self.noise.forward(frames)
     96 
     97         return frames, bg_frames

C:\Users\chris\.conda\envs\decode_env\lib\site-packages\decode\simulation\camera.py in forward(self, x, device)
     94         """Gamma for EM-Gain (EM-CCD cameras, not sCMOS)"""
     95         if self._em_gain is not None:
---> 96             camera = self.gain.forward(camera)
     97 
     98         """Gaussian for read-noise. Takes camera and adds zero centred gaussian noise."""

C:\Users\chris\.conda\envs\decode_env\lib\site-packages\decode\simulation\noise_distributions.py in forward(self, x)
     56 
     57     def forward(self, x):
---> 58         return torch.distributions.gamma.Gamma(x, 1 / self.scale).sample()
     59 
     60 

ZeroDivisionError: float division by zeroge value, simulation: {sim_frames.mean().round()}, real data: {data_frames.mean().round()}')

The "division by zero" suggested that DECODE had a problem with em_gain=0, so I replaced it by em_gain = 1 in the yaml file from SMAP:

Camera:
  baseline: 100.0
  e_per_adu: 0.48
  em_gain: 1.0
  px_size: [97.0, 97.0]
  read_sigma: 74.4
  spur_noise: 0.002
InOut: {calibration_file: 'D:/data SMLM/20210303 Calib 3D Karo/210209_2C_2Cameras_100nm/zStackCamera1_3dcal.mat',
  experiment_out: 'D:/data SMLM/20210303 Calib 3D Karo/210209_2C_2Cameras_100nm'}
Simulation:
  bg_uniform: [24.927, 38.2988]
  emitter_av: 25.0
  emitter_extent:
  - [-0.5, 39.5]
  - [-0.5, 39.5]
  - [-62.0, 507.0]
  intensity_mu_sig: [2284.5032, 456.90063]
  lifetime_avg: 1.0438

Now it runs up to the comparison between the simulated frames and the experimental frames. However, the simulated frame show only noise and no visible blinking event using the parameters from the SMAP file, with an intensity span from 0-300 on the simulation vs 0-80 on the experimental frame:

So there seem to be a problem here. I don't understand why I don't see any blinking event on the simulated frame with this yaml parameter file. If you want to reproduce using the yaml file above, I also attach the calibration file used:
http://www.neurocytolab.org/up/DECODE_sCMOS_Zcalib.mat

Thank you for your help!

PS.
A side issue:
The text just above the comparison step in the notebook says:

Looking at the mean brightness, we see that there is a large mismatch. The reason is that we simulate the photons that are emitted by the fluorophores, while the tiff file shows the photons that are recorded by the camera (after amplification). To enable direct comparison we convert the real data frames into photon numbers.
Be aware that camera.forward() and camera.forward() are not each others inverse. forward() performs noise sampling, while backward() simply rescales by the em_gain

Indeed, the mean intensity of the simulation is 44.0, while the average intensity of the experimental frames as reported in the notebook is 157.0 before the camera.backward rescaling and 27.0 after the rescaling (which is consistent with [(157-baseline)/(e_per_adu•em_gain)], with the values from the yaml file above (baseline 100, e_per_adu 0.48, em_gain 1). As an aside I think "rescaled by the em_gain" in the notebook test is imprecise as it's rescaled by (e_per_adu•em_gain) if I understand correctly what happens.

[Haydnspass]

Hey @cleterrier,

for no EM-Gain, it should be set to None not 0 or 1, but you have to explicitly set it to nothing, you must not leave it out.
Therefore this part should look like:

Camera:
  baseline: 100.0
  e_per_adu: 0.48
  em_gain: 
  px_size: [97.0, 97.0]
  read_sigma: 74.4
  spur_noise: 0.002

I'll have a look at your calibration file. If you want to debug you can overwrite your camera in the notebook with

# ...
camera = decode.simulation.camera.PerfectCamera(device='cpu')

which will only result in shot noise.

[cleterrier]

I added the perfect camera in the training notebook, just before:
# finally we derive some parameters automatically for easy use param = decode.utils.param_io.autoset_scaling(param)

and the result is the same pattern of pure shot noise. So I guess the Z calibration is the culprit. This is a spline calibration made by SMAP using a larger window (31x31 pixels) because there is a lot of astigmatism on this setup, so PF are extended. Here's the SMAP Matlab Figure corresponding to the calibration:
http://www.neurocytolab.org/DECODE_sCMOS_CalibFig.fig

and the corresponding Z stack tif (34 steps, 50 nm Z steps, 97 nm/pixel, 0.48 e-/ADU):
http://www.neurocytolab.org/DECODE_sCMOS_ZcalibStack.tif

[Haydnspass]

To me the PSF looks fine.
I sampled a few emitters and plotted a single frame. See the example below.
Good catch with the notebook. camera.forward sample from the noise distribution, camera.backward simply converts camera units back to photon units

[cleterrier]

Thanks @Haydnspass! Seems that the problem is not in the calibration then?

Can you try to execute the Training notebook with

• this calibration file (same as above)
• this yaml file from SMAP (with em_gain modified to none)
• these sample frames?

I'd be interested to see if you get something that looks like pure noise also... You can see the experimental data doesn't have a high SNR, maybe the SMAP parameters inferred make it even worse, and it just looks like noise?

Because with the above files, to me right now the comparison looks like this:

[Haydnspass]

Just tried, it is just the low photon count. As soon as I used a factor of 10 for intensity_mu_sig the emitters became visible to the eye (so no bug).
I checked the normalization of the PSF throughout the z range which looks fine as well. (see below, first normalization check, second one ROI without noise). Notice that the max. intensity is far less than the summed intensity so I guess intensity_mu is just off by a factor of 10-100.

Where did you get the intensity and bg values from? From the SMAP plugin or by guess?

[cleterrier]

I still don't understand why the simulation just shows noise visually in the Training notebook though - 2500 ± 500 photons per blinking events on a 30-50 background should be already very visible, and these the the values SMAP determines from the experimental sequence (where events are visible). I've thought about it while not sleeping last night - can't see any issue with gain value as sCMOS gain very straightforward (only one value for e/ADU, and em gain equivalent to 1) - also camera gain should not influence the simulation as intensity is in photons in the Figure shown. Background definition issue? Those mysterious camera parameters (74.4 and 0.002)? I mean as an SMLM microscopist, having to manually add a 10X factor to the number of photons (to 25,000) suggests there's a problem somewhere, as it's far from a typical Alexa Fluor 647 value.

One thing I can imagine is that as the astigmatism is large, the PSF is spread over a large number of pixels, so each pixel as a small fraction of the photon count (your PSF above at Z=0 has only 2.5% of the total photons - 0.025 to a normalized 1 total, would be even smaller away from the focal plane). But still the simulation should look like the experimental frames as it's derived from the SMAP estimates?

I also thought maybe the background was a problem. But making the background from 0-50 to 1-2 does not change the result - a shot noise images with a span of around 250 photons from the intensity scale bar.

What about camera baseline? Is it in photons or ADU for DECODE (as the simulation is in photons)? In SMAP it's in ADU (100 ADU)

[Haydnspass]

Hm a few thoughts

I still don't understand why the simulation just shows noise visually in the Training notebook though - 2500 ± 500 photons per blinking events on a 30-50 background should be already very visible, and these the the values SMAP determines from the experimental sequence (where events are visible).

I think the parameter extraction is the source of the problem, to be honest. Just by rough calculation: The frame on the right shows values of max ~250 (in photon units). Comparing this to the ROI I extracted from your PSF which has a max of 0.02 for 1 hypothetical photon this should (really) roughly correspond to around 10 000 photons [(250- bg:50)/0.02], given it was at the same position in z. Or have I made a mistake here?
With intensity_mu_sig ~2300+/-500 you would never end up with values like this since it's far beyond 3 sigmas off to the mean.

I've thought about it while not sleeping last night - can't see any issue with gain value as sCMOS gain very straightforward (only one value for e/ADU, and em gain equivalent to 1). Background definition issue? Those mysterious camera parameters (74.4 and 0.002)?

74.4 is gaussian noise, you can set it to 0 but I doubt it has a large effect. To simplify I would use PerfectCamera() until we know what the problem is here.

[Haydnspass]

Could you send me the localisations from SMAP extracted by the prefit? (You don't have to do it here publically, via Mail to [email protected] is okay) Then I can check the statistics

[cleterrier]

No problem to post everything here. I posted the files necessary to run the training above (3Dcal.mat, yaml out of SMAP, extracted images for comparison).

Here is the .sml.mat file from which the SMAP yaml file was generated (used fixed grouping at 97 nm to get the same yaml statistics as mine):
https://filesender.renater.fr/?s=download&token=ba9eec5e-8732-4466-acf5-d88cf7f88b48

Here is the tif sequence (768x768, 5K frames from a larger sequence) from which the localization .sml.mat file was calculated using SMAP (peak detection DoG 1.7, dynamic factor 1.3, window for detection 21 pixels, peak fitting Spline with 3dcal.mat above, window for fitting 21 pixels)
https://filesender.renater.fr/?s=download&token=1217739a-b528-4ced-942e-10f13532093d

Here is the cameras.mat file with the camera parameters used to process the sequence and obtain the sml.mat file (I find working with several cameras quite confusing, it keeps reverting to default, so I use a cameras.mat file with only one camera as default):
https://filesender.renater.fr/?s=download&token=64c4590d-b0b8-43fb-94a4-41241059192a
Should be 97 nm/pixel, 0.48 e-/ADU, no EM gain

(btw just to avoid misunderstanding, I use the edit function extensively in github to add precision and correct stuff on my messages so it's best to read here rather than using the email sent with the first "draft" version :)

[cleterrier]

OK I made some progress!
The camera.read_sigma was fixed at 74.4 from the SMAP yaml file. Looking at the spec sheet from my sCMOS camera, this is the read noise profile:

So the readout noise is at 1.65 e- with a sigma of around 0.5 e-
If I put camera.read_sigma at 0.5, all other parameters equal, I get something much closer between simulation and experimental frame (keeping the 2500±500 values for blinking events given by SMAP):

I guess the 74.4 value is after EM gain for an EMCCD? Not sure what are camera.read_sigma and camera.spur_noise exactly, but it looks like sigma should be much lower than 74.4 for an EMCCD and that the SMAP parameter evaluation should reflect this difference

[Haydnspass]

Yeah, I just had another look at the camera pipeline, the readout noise is added after em_gain in electron units.
Read sigma is the spread of the gaussian that is added (see below), therefore I think it must be 1.65 in your case. We do not account for the error of the error so to say.

I think you are still a bit low in photon counts

Spur_noise is spurious noise, but it won't make a difference.
You can have a look at the pipeline here:

DECODE/decode/simulation/camera.py

Lines 88 to 107 in 1aba0e4

	"""Clamp input to 0."""
	x = torch.clamp(x, 0.)
	"""Poisson for photon characteristics of emitter (plus autofluorescence etc."""
	camera = self.poisson.forward(x * self.qe + self.spur)
	"""Gamma for EM-Gain (EM-CCD cameras, not sCMOS)"""
	if self._em_gain is not None:
	camera = self.gain.forward(camera)
	"""Gaussian for read-noise. Takes camera and adds zero centred gaussian noise."""
	camera = self.read.forward(camera)
	"""Electrons per ADU, (floor function)"""
	camera /= self.e_per_adu
	camera = camera.floor()
	"""Add Manufacturer baseline. Make sure it's not below 0."""
	camera += self.baseline
	camera = torch.max(camera, torch.tensor([0.]).to(camera.device))

And that is how we convert from camera units back to photons

DECODE/decode/simulation/camera.py

Lines 131 to 137 in 1aba0e4

	out = (x - self.baseline) * self.e_per_adu
	if self._em_gain is not None:
	out /= self._em_gain
	out -= self.spur
	out /= self.qe
	return out

[cleterrier]

OK I had forgotten to comment out a perfect camera somewhere, so my backward conversion was wrong. Now it's perfect

[cleterrier]

This is using all the parameters from the SMAP yaml file, adjusting the em_gain (from 0 to no value), camera QE (from 1.0 to 0.8) and read_sigma (from 74.4 to 1.65).

[Haydnspass]

Great! :) So than that's solved, right?
Not much we can do about the camera can we? Or was it detected wrongly by SMAP?

[cleterrier]

Looks like it is! I think we learned that the camera noise value was important - could be good to reflect that in the docs. I'll check in SMAP where the 74.4 value comes from.

[Haydnspass]

@cleterrier Another thing that popped into my mind: Since you were talking about using SCMOS, have you turned off Mirror2D?
I have added a description of when and when not to use it here in the most recent Fitting notebook:
https://github.com/TuragaLab/DECODE/blob/master/decode/utils/examples/Fit.ipynb