Scaling of initial K0 guesses with TIFF files

[thomasmfish]

We are trying to compare cudasirecon output with that from softworx on some data we have already reconstructed with softworx. The cudasirecon is failing in a strange manner: we get a successful reconstruction but it appears to find stripes with exactly twice the expected line spacing. The the reconstructions have a small enhanced resolution but not the expected successful reconstruction.

We believe that the issue is that the initial k0 estimate is incorrectly calculated and the initial guess for the k0 position in pixels in Fourier space is incorrect. The angle is fine but the magnitude is half as big as expected.

I'm tagging @iandobbie as I'm working with him on this.

cudasirecon log (525nm)
softworx log

[tlambert03]

can you share the parameters you used as well? It's been a long time since I used softworx, and I can't remember off the top of my head, but it's conceivable that the spacings are expressed differently. I can tell you that on my OMX, for 488nm excitation, i use

k0angles=-0.804300,-1.8555,0.238800
ls=0.2035

(and in the log output for that wavelength, it settles on:

Optimum k0 angle=-0.803223, length=98.784332, spacing=0.414641 microns

so if you're starting with linespacing 0.39600 (which I pulled from your softworks log file), you might need to halve that.

maybe @linshaova, remembers whether there was a difference in the way GE expressed linespacing?

[thomasmfish]

Quick reply, thank you!

I've attached the config file (without the modified ls):
config525.txt

Why is the spacing 0.414641 if your input is 0.2035? We tried halving the line spacing and this caused it to fail (bad k0 fit) - I've just run it with a halved ls again, with the same file, to illustrate via the log below.

Halved line spacing log file (ls=0.197):
data525_recon_525.log

[tlambert03]

Why is the spacing 0.414641 if your input is 0.2035?

it is indeed confusing! 😂 I need to look back into the source code a bit to figure that seeming contradiction out.

The source code is here if you'd like to look yourself:

cudasirecon/src/cudaSirecon/gpuFunctionsImpl.cu

Lines 874 to 995 in 97af94d

	__host__ void fitk0andmodamps(std::vector<GPUBuffer>* bands,
	GPUBuffer* overlap0, GPUBuffer* overlap1, int nx, int ny, int nz,
	int norders, vector *k0, float dxy, float dz, std::vector<GPUBuffer>* otf,
	short wave, cuFloatComplex amps[], ReconParams * pParams)
	{
	int fitorder1 = 0;
	int fitorder2 = 0;
	if (nz > 1) {
	fitorder2 = 2;
	if (pParams->bBessel)
	fitorder2 = 1;
	}
	else {
	fitorder2 = 1;
	}
	float k0mag = sqrt(k0->x * k0->x + k0->y * k0->y); // in 1/um
	float k0angle = atan2(k0->y, k0->x);
	/* recalculate the overlap arrays at least this first time */
	int redoarrays = (pParams->recalcarrays >= 1);
	float x2 = k0angle;
	cuFloatComplex modamp;
	float amp2 = getmodamp(k0angle, k0mag, bands, overlap0, overlap1, nx, ny, nz,
	fitorder1, fitorder2, dxy, dz, otf, wave, &modamp, redoarrays, pParams, 0);
	/* recalculate the overlap arrays every time only if recalcarrays >= 3*/
	redoarrays = (pParams->recalcarrays >= 3);
	float deltaangle = 0.001;
	float deltamag = 0.1 / (std::max(nx, ny) * dxy); // in 1/um
	float angle = k0angle + deltaangle;
	float x3 = angle;
	float amp3 = getmodamp(angle, k0mag, bands, overlap0, overlap1, nx, ny, nz,
	fitorder1, fitorder2, dxy, dz, otf, wave, &modamp, redoarrays, pParams, 0);
	float amp1;
	float x1 = 0.0;
	float a;
	if (amp3 > amp2) {
	while(amp3 > amp2) {
	amp1 = amp2;
	x1 = x2;
	amp2 = amp3;
	x2 = x3;
	angle += deltaangle;
	x3 = angle;
	amp3 = getmodamp(angle, k0mag, bands, overlap0, overlap1, nx, ny, nz,
	fitorder1, fitorder2, dxy, dz, otf, wave, &modamp, redoarrays, pParams, 0);
	}
	} else {
	angle = k0angle;
	a = amp3;
	amp3 = amp2;
	amp2 = a;
	a = x3;
	x3 = x2;
	x2 = a;
	while (amp3 > amp2) {
	amp1 = amp2;
	x1 = x2;
	amp2 = amp3;
	x2 = x3;
	angle -= deltaangle;
	x3 = angle;
	amp3 = getmodamp(angle, k0mag, bands, overlap0, overlap1, nx, ny, nz,
	fitorder1, fitorder2, dxy, dz, otf, wave, &modamp, redoarrays, pParams, 0);
	}
	} /* the maximum of modamp(x) is now between x1 and x3 */
	angle = fitxyparabola(x1, amp1, x2, amp2, x3, amp3); /* this should be a good angle. */
	/***** now search for optimum magnitude, at this angle *****/
	x2 = k0mag;
	amp2 = getmodamp(angle, k0mag, bands, overlap0, overlap1, nx, ny, nz,
	fitorder1, fitorder2, dxy, dz, otf, wave, &modamp, redoarrays, pParams, 0);
	float mag = k0mag + deltamag;
	x3 = mag;
	amp3 = getmodamp(angle, mag, bands, overlap0, overlap1, nx, ny, nz,
	fitorder1, fitorder2, dxy, dz, otf, wave, &modamp, redoarrays, pParams, 0);
	if (amp3 > amp2) {
	while (amp3 > amp2) {
	amp1 = amp2;
	x1 = x2;
	amp2 = amp3;
	x2 = x3;
	mag += deltamag;
	x3 = mag;
	amp3 = getmodamp(angle, mag, bands, overlap0, overlap1, nx, ny, nz,
	fitorder1, fitorder2, dxy, dz, otf, wave, &modamp, redoarrays, pParams, 0);
	}
	} else {
	mag = k0mag;
	a = amp3;
	amp3 = amp2;
	amp2 = a;
	a = x3;
	x3 = x2;
	x2 = a;
	while (amp3 > amp2) {
	amp1 = amp2;
	x1 = x2;
	amp2 = amp3;
	x2 = x3;
	mag -= deltamag;
	x3 = mag;
	amp3 = getmodamp(angle, mag, bands, overlap0, overlap1, nx, ny, nz,
	fitorder1, fitorder2, dxy, dz, otf, wave, &modamp, redoarrays, pParams, 0);
	}
	} /* the maximum of modamp(x) is now between x1 and x3 */
	mag = fitxyparabola(x1, amp1, x2, amp2, x3, amp3); /* this should be a good magnitude. */
	/* if we were perfectionist we'd iterate for angle again now */
	printf("Optimum modulation amplitude:\n");
	redoarrays = (pParams->recalcarrays>=2); /* recalculate the d_overlap arrays for optimum modamp fit */
	amp3 = getmodamp(angle, mag, bands, overlap0, overlap1, nx, ny, nz,
	fitorder1, fitorder2, dxy, dz, otf, wave, &modamp, redoarrays, pParams, 1);
	/* one last time, to find the modamp at the optimum k0*/
	printf("Optimum k0 angle=%f, length=%f, spacing=%f um\n", angle, mag, 1.0 / mag);

I'll try to find out something more concrete for you, and maybe @linshaova will also have thoughts

[linshaova]

Hi everyone,

Sorry for the confusion. Line spacing output from the code should be roughly twice the number you input with --ls flag in 3D SIM. The output line space corresponds to the lowest order lateral SIM pattern's line spacing, which in 3D SIM is double the finer pattern. This was due to our effort to unify the code so that line spacing output (not input) has the same definition for all modes of SIM: 3D SIM, 2D SIM, and nonlinear SIM.

lin

[tlambert03]

thanks @linshaova.

@thomasmfish, I just double checked my softworx settings, and they are indeed very close to my cudasirecon settings (for dv files)... on softworx as well, I have a linespacing of ~208nm for the 488x/525m channel. So, I'm a bit surprised that your line spacing would be upwards of 400nm even on softworx. Are you doing anything unusual that would cause your lines to be spaced so far apart (like, are you using a low NA objective or doing something different?).

one thing I'm remember now... if I recall correctly, dv files have an inverted Y axis compared to tiff files, and if you look at the two config files in our test data here, the angles are indeed different:

cudasirecon/test_data/config

Lines 1 to 5 in 97af94d

	nimm=1.515
	background=0
	wiener=0.001
	k0angles=-0.804300,-1.8555,0.238800
	ls=0.2035

cudasirecon/test_data/config-tiff

Lines 1 to 6 in 97af94d

	nimm=1.515
	fastSI=1
	background=0
	wiener=0.001
	k0angles=0.804300,1.8555,-0.238800
	ls=0.2035

so, perhaps, drop your linespacing back down to the one closer to 200nm, and flip the sign on all of your angles

[thomasmfish]

@linshaova Thanks for your input. I believe we were expecting the 2nd order line spacing to be found and it is only finding the 1st - @iandobbie has a much deeper understanding than I do, so hopefully he can confirm whether I'm using the correct terms here. We were expecting double the magnitude for the k0s and were confused by why the line spacing wasn't consistent, so we thought there might be some scaling issue somewhere (perhaps just for tiffs - hence the name of the issue).

@tlambert03 We are using an NA of 0.9, and I've made sure the line spacings correctly match softworx for each wavelength, so I believe they should be correct.
The k0 signs are a good point - I hadn't considered that. It seemed to be getting good matches for the k0 angles but I don't believe I have tried inverting them. I am surprised it managed to fit so closely for such incorrect k0s if that is the case though. I'll have to try tomorrow and get back to you.

[tlambert03]

I am surprised it managed to fit so closely for such incorrect k0s if that is the case though.

yeah, i agree, that would be the surprising thing, but I've definitely seen it "settle" into incorrect local minima before when it was way off.

I believe we were expecting the 2nd order line spacing to be found and it is only finding the 1st

in any case, the number that the software is expecting you to enter is definitely that outer order, so, if you are using a 0.9 NA lens then your spacing in the upper 300nm range does sound about right. My guess/hope is that if you stick with that spacing, and invert the angles, it should be fine :) let us know how it goes

[thomasmfish]

Inverting the angles causes the reconstructions to fail, and I realised that we did check the magnitude against an FFT and found it matched the 1st order magnitude. At this point, I think the issue may be that it is failing to find the 2nd order, which is why it is blurrier than expected. Or it's a scaling factor. I'll spend some time looking through in more detail today.

[tlambert03]

Could you perhaps share the data?

One more question: since you're able to reconstruct this in softworx, you must be using a dv file there right? (I don't recall them supporting tiff?). Just want to know whether this is a parameter issue, or a file type issue

[iandobbie]

I think something in the code must be flipping the angles as the fits succeed and produce angles very similar to softworx and output spacings close to exactly twice the softworx output. The reconstructions we are getting end pretty sharply at about 300 nm rather than the 200 nm which made me think we were getting about 2^0.5 improvement from the coarse stripes rather than the 2x expected from a full 3d SIM reconstruction.

I guess out issue must be something else, maybe the OTFs? there is a bit of variability but the fitting amplitudes look similar to the softworx results so I am surprised about the pretty dramatic reduction of the reconstruction resolution.

Another avenue I have investigated but not got to the bottom of is possible differing scale of the Wienner parameter so we are effectively filter much more aggressively with the cuda code compared to the softworx reconstruction.

Thanks for your help and pointers so far. We are getting towards a working pipeline with much faster reconstructions than previously.

[iandobbie]

Could you perhaps share the data?

Tom will have to answer the data sharing issue, not mine to share I'm afraid

One more question: since you're able to reconstruct this in softworx, you must be using a dv file there right? (I don't recall them supporting tiff?). Just want to know whether this is a parameter issue, or a file type issue

The original data are .dv files. Tom has code that is taking the .dv files and dumping them into the recon process based on the python code in the repository but with adaptions to suit their workflow. The issue may well be file type related and I tried to build the cuda code to work with .dv files but struggled with cuda library versions etc... I am traveling for a couple of weeks and only have a mac laptop with no nvidia card so can work on this until I return.

[tlambert03]

I think something in the code must be flipping the angles as the fits succeed and produce angles very similar to softworx and output spacings close to exactly twice the softworx output.

The reason I don't think this is the case is "generally" is that I'm using very similar parameters (with dv files) in both softWoRx and cudasirecon. (Granted it's been a long time since I used softWoRx, but when I checked yesterday, the params I was using were the same as what I used to use). And when I did my original pipeline switch as you are doing, I eventually found extremely similar reconstruction results once I got the parameters right.

However, I still almost exclusively use dv files and not tiff.

That's why I'm currently less convinced by the theory that the code is wrong, and more curious about what might be different about your use case or file (eg tiff vs dv). I've never tried an air objective with low NA though, so it's possible that there's an error along those lines that I hadn't run into

Might need to see the data to be of much more help

[thomasmfish]

I have been given permission to share the data, though I'm not sure about posting it publicly. How should I get it to you?

As Ian said, we are using DV files in softworx. For cudasirecon, I am splitting the DV files by channel, saving to TIFFs and config files (the config containing settings from the DV metadata), then passing to pycudasirecon via numpy array.

[tlambert03]

You can email a link? Talley.lambert@gmail.com

I'll mention the following codebase as well:
https://github.com/tlambert03/otfsearch

The code there is absolutely hideous 🤢 😅 and 10 years old (using python 2). I wrote it before I created the pycudasirecon wrapper and never got around to updating... but that's what all my users use on our OMX to split channels/try a bunch of OTFs, pick the best one, reconstruct and merge them back together, etc... it stays with dv/mrc files all the way through (but conceptually, it should be fine with tiff as well, we just didn't have tiff support when I wrote that). Might find something to compare to in there amidst all the flotsam :)

[tlambert03]

Also inasmuch as you're ultimately using pycudasirecon, you might consider reading and splitting your dv files directly into numpy buffers (without writing to tiff intermediate) using mrc: https://github.com/tlambert03/mrc