Check & revise LN2_DEVEIN linear mode

[ofgulban]

A reminder for myself to revise LN2_DEVEIN 's -linear mode. The output value range seems problematic due to the scaling.

[ppxdm4]

Hi @ofgulban, I've been looking at this option of LN2_DEVEIN recently and think I'm having the same issue. The 'deveined' profile looks different when you use layer depths 0->1 (metric option) compared to 1->6 (for nr_layers=6). The non-metric case is corrected the profile such that it has a very steep negative gradient. In the code snippet below (lines 243-259) I believe l = 1 - *(nii_layer_data + i) * layer_max; should actually be l = 1 - *(nii_layer_data + i) / layer_max;

if (mode_linear) {  // Handle linear case straightforwardly
        for (int t=0; t<size_t; ++t) {
            for (int i=0; i<nr_voxels; ++i) {
                int j = t * nr_voxels + i;
                if (*(nii_layer_data + i) > 0) {
                    float l;
                    if (layer_max <= 1.0) {  // Indicates metric file is used
                        l = 1 - *(nii_layer_data + i);
                    } else {
                        l = 1 - *(nii_layer_data + i) * layer_max;
                    }

                    *(nii_output_data + j) = *(nii_input_data + j) * l;
                }
            }
        }
    }

Apologies if this is incorrect or you have already addressed the issue. It seems to have worked with the simple cases I have been testing on and makes sense when I think about how variable l is used.

[ofgulban]

Hi @ppxdm4 , Thanks a lot for bringing this to our attention and your informative comment. We ended up not using LN2_DEVEIN over the last year and were postponing looking into this issue. Together with your comment, now we have a reason to address this issue (tagging @layerfMRI as this might be of interest to him too).

[layerfMRI]

Thanks @ppxdm4, your fix is now included in the devel branch. I added another shift of half a layer, to account for the fact that there is no zeroth layer (the average centroid of layer 1 is at 0.5). And it looks as expected in a quick test.

It is only now that I appreciate what @ofgulban did in 1d5a3d95eaf920196674f2c5dce6a5f129f95014 (related to issue #28). It's not a scaling with the (multiplicative) reciprocity of the layers, but a multiplication with the (additive) inverse of the depth. So my earlier comment during the second coffee break is no longer valid. I think this solves my issue about normalizing it to layer 1. I think this issue can be closed at the next merge to master

[ppxdm4]

Hi @ofgulban, @layerfMRI, thank you both for these updates. I have pulled the latest devel branch and will be implementing the new linear deveining in the next few days. I have a quick question after reading issue #28: Could you clarify - when using LN2_DEVEIN with the -linear flag, if I provide a columns file will the corrections be performed on a column by column basis?

[layerfMRI]

Hi @ppxdm4, Yes, I believe that the linear case does not use the column file, if you provide it or not. The linear scaling happens on a voxel by voxel level. This is different to the deconvolution case. Which works on a column by column basis. This was different before @ofgulban revised it in 1d5a3d9. I am not sure, why this was changed, but I think the results should be comparable for thin columns. The new way is more minimalistic (more safe). The new way does not average within the columns. Only this way it is possible to straightforwardly use the metric file too. The new way of ignoring the columns for the linear case does not have the SNR benefit of local averaging within the layer-column intersection. And thus, it is not easily comparable with the deconvolution case.

[ofgulban]

Looking at the previous comments on #28, I think we have changed the columns file behavior due to a few discrepancies between the code and the blog post. Looking at the code, yes the columns file seems to be not used with linear option right now (v2.0.0).

However, this is not a strong limitation. We can code in so that linear option uses columns to first averages the voxels within each column, does linear deveining, and writes back into the columns. If @ppxdm4 want to use linear deveining in such a way, I am happy to look into this in the upcoming days :)