ENH dwidenoise: Overcomplete local PCA

[Lestropie]

https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0073021

Currently, in dwidenoise, the reconstructed signal for a given voxel is determined exclusively by the reduced-rank PCA decomposition from the sliding spatial window kernel when centred at that voxel.

What can be done instead is that the output intensities for a given voxel can be a weighted average of the rank-reduced decompositions of multiple positions of the sliding spatial window kernel; ie. not exclusively when the sliding window is centred at that same voxel. In the article cited above these weights are determined based on the number of components of each decomposition; I could imagine hypothetically that spatial proximity of the output reconstructed voxel to the centre of the kernel could also be used. This might make denoised outputs less sensitive to threshold effects on preservation / rejection of particular components.

Would need to either:

• Mutex across threads to control updating of output image intensities;
• Each thread would need a scratch buffer of the same size as the output series, updating with those patches processed by that thread, with the output image finally being formed from the summation of such.

The former is probably safer as RAM requirements could get quite large with the latter.

[Lestropie]

Given I've gotten myself hooked on denoising improvements over the last week, I think I'd like to take this on.

I am however interested to know if anybody has thoughts on what the recombination should look like. For each output image voxel, there is a set of denoising patches of which that voxel is a member. For each of those patches, there will be both an estimated signal rank, and a distance between the centre of the patch and the voxel. So one needs some algebraic expression that aggregates the reconstructed output intensities for that voxel from the different patches:

• In the original article linked above these recombination weights are based on the ranks of the patches of which that voxel is a member (see Equation 3 there).
• In the tensor MPPCA denoising article they mostly advocate for equal weighting of all patches of which the voxel is a member, though their discussion was not conclusive.
• In Lucilo's work, a Gaussian kernel was used for simplicity, though I don't see any information on the FWHM of the kernel used.

I suspect there will need to be support for multiple such expressions.

Also, in the scenario where no patch aggregation is performed, ideally the code should persist with the current optimisations where only the PCA projection of the signal in the voxel at the centre of the patch is computed, whereas for overcomplete variants the PCA projection for the whole patch will need to be computed.