Bug/groupwise registration

[bpinsard]

It seems that the way time interpolation is carried is wrong for interleaved at least.
So if I got it right, what we want from this function (mainly for scanner_time in Image4D) is an offset in time coordinates for sampling using 4D cubic splines.
here is an example of what it does:

>tr=3
>n_slices=48
>slice_order=range(0,n_slices,2)+range(1,n_slices,2)
>interp = gr.interp_slice_order(np.arange(0,n_slices),slice_order)
>interp
array([ 0,  2,  4,  6,  8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32,
       34, 36, 38, 40, 42, 44, 46,  1,  3,  5,  7,  9, 11, 13, 15, 17, 19,
       21, 23, 25, 27, 29, 31, 33, 35, 37, 39, 41, 43, 45, 47])
>interp*float(n_slices)
array([ 0.        ,  0.04166667,  0.08333333,  0.125     ,  0.16666667,
        0.20833333,  0.25      ,  0.29166667,  0.33333333,  0.375     ,
        0.41666667,  0.45833333,  0.5       ,  0.54166667,  0.58333333,
        0.625     ,  0.66666667,  0.70833333,  0.75      ,  0.79166667,
        0.83333333,  0.875     ,  0.91666667,  0.95833333,  0.02083333,
        0.0625    ,  0.10416667,  0.14583333,  0.1875    ,  0.22916667,
        0.27083333,  0.3125    ,  0.35416667,  0.39583333,  0.4375    ,
        0.47916667,  0.52083333,  0.5625    ,  0.60416667,  0.64583333,
        0.6875    ,  0.72916667,  0.77083333,  0.8125    ,  0.85416667,
        0.89583333,  0.9375    ,  0.97916667])

which means that the second slice for example, acquired at half tr is applied a very small offset.
with this fix it gives:

>interp = gr.interp_slice_order(np.arange(0,n_slices),slice_order)
>interp
array([ 0, 24,  1, 25,  2, 26,  3, 27,  4, 28,  5, 29,  6, 30,  7, 31,  8,
       32,  9, 33, 10, 34, 11, 35, 12, 36, 13, 37, 14, 38, 15, 39, 16, 40,
       17, 41, 18, 42, 19, 43, 20, 44, 21, 45, 22, 46, 23, 47])
>interp/float(n_slices)
Out[207]: 
array([ 0.        ,  0.5       ,  0.02083333,  0.52083333,  0.04166667,
        0.54166667,  0.0625    ,  0.5625    ,  0.08333333,  0.58333333,
        0.10416667,  0.60416667,  0.125     ,  0.625     ,  0.14583333,
        0.64583333,  0.16666667,  0.66666667,  0.1875    ,  0.6875    ,
        0.20833333,  0.70833333,  0.22916667,  0.72916667,  0.25      ,
        0.75      ,  0.27083333,  0.77083333,  0.29166667,  0.79166667,
        0.3125    ,  0.8125    ,  0.33333333,  0.83333333,  0.35416667,
        0.85416667,  0.375     ,  0.875     ,  0.39583333,  0.89583333,
        0.41666667,  0.91666667,  0.4375    ,  0.9375    ,  0.45833333,
        0.95833333,  0.47916667,  0.97916667])

which seems more correct.

also interpolation between 2 slices having half tr offset

>interp = gr.interp_slice_order(np.arange(0,1.1,.1),slice_order)
>interp/float(n_slices)
array([ 0.  ,  0.05,  0.1 ,  0.15,  0.2 ,  0.25,  0.3 ,  0.35,  0.4 ,
        0.45,  0.5 ])

is this correct?

[alexis-roche]

Thanks for checking that up. You seem to be misunderstanding the meaning of the slice_order argument. In this example, the second physical slice (from bottom to top of the head) is acquired third in time (slice_order[1]==2), hence much earlier than at half tr. The interleaved acquisition scheme you seem to have in mind is not the one encoded by the particular slice_order you defined.

The purpose of the interp_slice_order function is to interpolate the acquisition time (in slices counts) of fictional slices based on their respective real adjacent slices. So, if the Z argument consists of integer values smaller than the number of slices, meaning that all our fictional slices happen to be real, then the function should return the corresponding slice order values without interpolation. For instance, for Z==1 we should get 2 here, not 24.

[bpinsard]

ok so the interleaved acquisition scheme that you defined in Image4D would acquire in the following order:
-slice bottom
-slice middle of the volume
-slice bottom+1
-slice middle of the volume+1
etc...
which is not something really common, or I am wrong again?
for me the interleaved scheme would be:
slice 0
slice 2
....
slice 1
slice 3

[bpinsard]

In fact this mean there would be a difference between the spm type of slicetiming definition that is listing the indices of the slices being acquired in temporal order and yours which seems to be defining the slice_order as the list of temporal position sorted by slice order.
Do you see what I mean?

[alexis-roche]

I do see what you mean and I think this discrepancy between SPM and Realign4D regarding slice timing definition was a source of confusion in a previous discussion:

nipy#134

My initial default interleaved slice order was indeed what you have in mind (say slice_order = [0, 4, 1, 5, 2, 6, 3] meaning [0, 2, 4, 6, 1, 3, 5] in the SPM convention). But then I changed it after being told by @matthew-brett that this was not the common interleaved acquisition scheme people would think of. From your comment, I am getting to realize that I got confused in that discussion and that Matthew was probably thinking of slice_order in the same way as you and most SPM users would. So my initial choice was actually correct in terms of implementing the "standard" interleaved scheme, and I changed it to something non-standard.

But most likely, I should re-define the slice_order argument to make it consistent with SPM and avoid that sort of confusion.

Very good point, thanks!

[bpinsard]

ok, I think due to the fact that the default interleaved scheme of Image4d is not the usual expected for interleaved we can expect that people have been using it thinking of regular interleaved, maybe we should broadcast this information on the nipy list ?

[alexis-roche]

Please don't hesitate to start a thread on this. Otherwise I'll do it ASAP.

Alexis

On 20 mars 2013, at 15:49, bpinsard [email protected] wrote:

ok, I think due to the fact that the default interleaved scheme of Image4d is not the usual expected for interleaved we can expect that people have been using it thinking of regular interleaved, maybe we should broadcast this information on the nipy list ?

—
Reply to this email directly or view it on GitHub.