Idea: use as a simulator for user training ?

[NicoKiaru]

First of all, super nice project!

Quick question, have you by any chance considered this repo as a way to simulate a 'real' live confocal acquisition ? My use case would be to train people without a microscope. Controls for focus, scan speed, pinhole, etc would appear and update live, like if you were sitting in front of a microscope.

Thanks!

[tlambert03]

hey @NicoKiaru

absolutely. However... i don't think the speed would (currently) be satisfying enough for the use case of having a GUI with the usual acquisition settings, with live updates as you change them. If we were strictly discussing (slow) point scanners, and if we strictly limited the output to a single plane (rather than the full 3D stack), then we could probably make some single-plane-specific speed improvements here, and possibly match the ~1-2 seconds of a real point scanner, assuming you were using one of the GPU backends.

[NicoKiaru]

Oh whoa, I wasn't aware that this was so compute-intensive. Thanks for the fast answer. I hope I find some time to play with microsim.

[tlambert03]

the compute intensiveness ultimately comes down to (and can be balanced with) simulation accuracy. The biggest cost is 3D convolutions of large volumes with large PSF kernels.

In a "worst-case" scenario (which we'll relax in a moment), imagine we're starting with the segmentations from the Cosem/OpenOrganelle datasets as ground truth. Those are sampled at 4nm voxels... so if you wanted to make a final image with a field of view that encompasses the entire cell/dataset, you're starting with ~1000 x 12,000 x 12,000 voxels. And, then if you want to make a "realistic" PSF that accurately models the spread of out-of-focus light into other planes. You need a PSF that is also pretty big (up to 1000s+ pixels per side) if you don't want to have "sharp edges" at certain points in the resulting convolution (provided sufficient SNR). So, that convolution alone would be crippling. So, one needs to make concessions, either to the field of view (in X,Y,or Z) or the "realism" of the ground truth, by downsampling to ground truth voxel sizes of 32nm or even more (and fluorophores are clearly not distributed as blocky cubes of 32nm^2), or by cutting down the PSF size and living with edge artifacts, or by "hiding" those artifacts by adding a large amount of detector noise.

Point is, there are many decisions to be made... microsim does let you tune all of those things relatively easily... so I do think it would be pretty fun and possible to make a GUI that restricts you to a relatively responsive range of parameters. And I also have little doubt that, for the purposes of education, many of the compromises made to accuracy there would be largely negligible.

So, long story short, I think it's a great idea 😂 ... we'll just see how far we need to go to decrease the complexity/realism to make it responsive

---

edit: even if one is ultimately after a 2D output image... at least a small 3D volume as input must be provided if we are to accurately show the effect of, say, pinhole size on rejecting out of focus light... since we need to have provided that out-of-focus light emitter to begin with

[NicoKiaru]

Nice! Thanks for the explanation. Indeed, I see how giving away precision can lead to strong speed-ups. Are ground-truth dataset multi-resolution already ? Or does one need to compute downscaled blocks ?

[NicoKiaru]

Ah, cosem have multiple resolution, great!

[tlambert03]

yeah... but it's actually much more complicated than that when it comes to organelle segmentations (explanation below if you're curious). The good new is that microsim handles all that for you. All you need to do is specify the scale of the ground-truth space as something like 0.032 and microsim takes care of the rest. See, for example, the cosem example here... where we start with 64nm ground truth in the interest of showing something quickly :)

---

nitty gritty details:

• OpenOrganelles primary goal was instance segmentation, not just semantic pixel classification (which is lovely)... and therefore the values of their label maps are not just 0 or 1, but rather 0-N where N is the number of objects.
• In their downscaled versions, they chose to go with the mode value of the pixels in the parent scale... (rather than the mean or sum). This non-linearity was introduced so that you could still load the downscaled version in something like neuroglancer and see something like "this pixel was mostly made up of mitochondria number N". However... for our purposes, that's not useful, because we want the lower resolutions to represent the density/abundance of molecules in that voxel. So, in microsim, when you use a CosemDataset label with a ground-truth scale of less than 0.004 (raw data), it will download the scale0 data, rebin it using summing, and cache that locally.

this was mostly done in #50 ... though there's not any more detail in the PR description.

[tlambert03]

if you want to play directly with the cosem data... use microsim.cosem.CosemDataset

In [4]: from microsim.cosem import CosemDataset

In [5]: ds = CosemDataset.fetch('jrc_hela-3')

In [7]: print([img.name for img in ds.segmentation_layers])
[
    'cent-dapp_seg',
    'cent_seg',
    'chrom_seg',
    'ecs_seg',
    'endo-mem_seg',
    'endo_seg',
    'er-mem_seg',
    'er_seg',
    'eres_seg',
    'golgi-mem_seg',
    'golgi_seg',
    'gt',
    'ld-mem_seg',
    'ld_seg',
    'lyso-mem_seg',
    'lyso_seg',
    'mito-mem_seg',
    'mito_seg',
    'mt-out_seg',
    'ne-mem_seg',
    'ne_seg',
    'nhchrom_seg',
    'np_seg',
    'nucleolus_seg',
    'nucleus_seg',
    'pm_seg',
    'ribo_seg',
    'vesicle-mem_seg',
    'vesicle_seg'
]


In [9]: ds.image(name='golgi-mem_seg')
Out[9]:

CosemImage(
    name='golgi-mem_seg',
    description='Golgi membrane',
    url='s3://janelia-cosem-datasets/jrc_hela-3/jrc_hela-3.n5/labels/golgi-mem_seg',
    format='n5',
    grid_scale=[3.24, 4.0, 4.0],
    grid_translation=[0.0, 0.0, 0.0],
    grid_dims=['z', 'y', 'x'],
    grid_units=['nm', 'nm', 'nm'],
    sample_type='label',
    content_type='segmentation',
    created_at=datetime.datetime(2022, 10, 28, 3, 14, 31, 836345),
    dataset_name='jrc_hela-3',
    stage='prod',
    scales=['s0', 's1', 's2', 's3', 's4']
)

In [11]: ds.image(name='golgi-mem_seg').read(level=4)  # jump straight to resolution scale 4
[12:34:47:284393] INFO     Skipping level 1 because it already exists.               _tstore.py:242
Processing /Users/talley/Library/Caches/microsim/janelia-cosem-datasets/jrc_hela-3/jrc_hela-3.n5/labels/golgi-mem_seg/sum2
100%|█████████████████████████████████████████████████████████| 1200/1200 [00:01<00:00, 802.51it/s]
Processing /Users/talley/Library/Caches/microsim/janelia-cosem-datasets/jrc_hela-3/jrc_hela-3.n5/labels/golgi-mem_seg/sum3
100%|███████████████████████████████████████████████████████████| 156/156 [00:00<00:00, 935.23it/s]
Processing /Users/talley/Library/Caches/microsim/janelia-cosem-datasets/jrc_hela-3/jrc_hela-3.n5/labels/golgi-mem_seg/sum4
100%|████████████████████████████████████████████████████████████| 42/42 [00:00<00:00, 2125.95it/s]
Out[11]:

TensorStore({
  'context': {
    'cache_pool': {'total_bytes_limit': 4000000000},
    'data_copy_concurrency': {},
    'file_io_concurrency': {},
    'file_io_locking': {},
    'file_io_memmap': False,
    'file_io_sync': True,
  },
  'driver': 'n5',
  'dtype': 'uint16',
  'kvstore': {
    'driver': 'file',
    'path': '/Users/talley/Library/Caches/microsim/janelia-cosem-datasets/jrc_hela-3/jrc_hela-3.n5/labels/golgi-mem_seg/sum4/',
  },
  'metadata': {
    'blockSize': [256, 256, 256],
    'compression': {'level': -1, 'type': 'gzip', 'useZlib': False},
    'dataType': 'uint16',
    'dimensions': [775, 62, 750],
    'pixelResolution': {
      'dimensions': [64.0, 64.0, 51.84000000000002],
      'unit': 'nm',
    },
    'transform': {
      'axes': ['z', 'y', 'x'],
      'order': 'C',
      'scale': [51.84000000000002, 64.0, 64.0],
      'translate': [24.30000000000001, 30.0, 30.0],
      'units': ['nm', 'nm', 'nm'],
    },
  },
  'transform': {
    'input_exclusive_max': [[62], [750], [775]],
    'input_inclusive_min': [0, 0, 0],
    'input_labels': ['y', 'x', 'z'],
    'output': [
      {'input_dimension': 2},
      {'input_dimension': 0},
      {'input_dimension': 1},
    ],
  },
})

[NicoKiaru]

Makes total sense, I see how a proper downscaling for instance segmentation requires a complicated representation, like a dict {id->sum_of_voxels}, very annoying to represent. I guess it's not far from what you did in the repo.