Well, that’s the thing, we may be talking about two different thing: the original koay paper mentions a few equations, which can be used to correct the bias for both the signal mena and the noise std using as a reference the magnitude SNR for each voxel. To do it properly as in the mppca mrm from 2015, you need, once you have split the eigen values, to use the lower half as the 3D noise map and apply it for correcting the signal on the rest of the denoised 4D maps. Doing so, you end u with a voxelwise 4D correction factor, which you can apply on both the 4D data and the 3D noise maps, although it would give you a stack of 4D noise maps after correction.
I guess I can live with that, although then in the end it would be more sensible from the physics point of view to have a single 3D noise map in the end, since it would mean each voxel through time is subject to the same contamination effet, instead of reducing back the 4D map with say the median value (although I would expect it to be pretty close in most cases), so maybe it’s a good enough workaround.
Can’t find the bessel function (they have the spherical version in eigen though, not sure how it applies versus the regular one) in those special functions’ libs, but anyway, I’ll stay with my last week quick R core wrapper, it’s overkill for now but it gets stuff done, thanks for the tip.