searchlight did not run through the edge of input data

[yushiangsu]

I encountered a problem that the brainiak.searchlight always return the output have None at the edge of data. This is problematic because I did searchlight on native EPI images with minimal preprocessing. Some participants have bigger brain that cover the top and bottom edge of native EPI image. These None were produced at top and bottom edge of brain, but these voxels should have eligible searchlight results.

I'm not sure if this is a bug. And I have not seen any warning about this in the document. I would like to post the issue here, hoping this could be fixed, or could be added as a warning in the document to inform people who are going to use searchlight module.

Here is a minimal duplicable example to present the issue:
I first generate a toy data and specify mask as all ones (include all data in searchlight).

# import package
import numpy as np
from brainiak.searchlight.searchlight import Searchlight
# making a toy data
dsize = 10
t = 20
data = np.zeros((dsize,dsize,dsize,t))
for i in range(dsize // 2):
    data[i:(dsize-i), i:(dsize-i), i:(dsize-i), :] += 1
mask = np.ones((dsize,dsize,dsize), dtype = int)
# print the toy data
print(f"data shape: {data.shape}")
print(data[dsize//2, :, :, 0])
print(f"mask shape: {mask.shape}")
print(mask[dsize//2, :, :])

data shape: (10, 10, 10, 20)
[[1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 2. 2. 2. 2. 2. 2. 2. 2. 1.]
 [1. 2. 3. 3. 3. 3. 3. 3. 2. 1.]
 [1. 2. 3. 4. 4. 4. 4. 3. 2. 1.]
 [1. 2. 3. 4. 5. 5. 4. 3. 2. 1.]
 [1. 2. 3. 4. 5. 5. 4. 3. 2. 1.]
 [1. 2. 3. 4. 4. 4. 4. 3. 2. 1.]
 [1. 2. 3. 3. 3. 3. 3. 3. 2. 1.]
 [1. 2. 2. 2. 2. 2. 2. 2. 2. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]]
mask shape: (10, 10, 10)
[[1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]
 [1 1 1 1 1 1 1 1 1 1]]

Creating the simplest function that apply on each searchlight

# this will return the centroid of the searchlight
def sl_func(sldata, slmask, myrad, bcvar):
    return sldata[0][myrad, myrad, myrad, 0]

Implementing searchlight with radius = 3 and print the output

rad = 3 # specify searchlight radius
sl = Searchlight(sl_rad = rad, max_blk_edge=10)
sl.distribute([data], mask)
sl.broadcast(0)
sl_results = sl.run_searchlight(sl_func, pool_size = 1)
print(f"output shape: {sl_results.shape}")
print(sl_results[dsize//2, :, :])

output shape: (10, 10, 10)
[[None None None None None None None None None None]
 [None None None None None None None None None None]
 [None None None None None None None None None None]
 [None None None 4.0 4.0 4.0 4.0 None None None]
 [None None None 4.0 5.0 5.0 4.0 None None None]
 [None None None 4.0 5.0 5.0 4.0 None None None]
 [None None None 4.0 4.0 4.0 4.0 None None None]
 [None None None None None None None None None None]
 [None None None None None None None None None None]
 [None None None None None None None None None None]]

I tried with another radius = 2 and print the output

output shape: (10, 10, 10)
[[None None None None None None None None None None]
 [None None None None None None None None None None]
 [None None 3.0 3.0 3.0 3.0 3.0 3.0 None None]
 [None None 3.0 4.0 4.0 4.0 4.0 3.0 None None]
 [None None 3.0 4.0 5.0 5.0 4.0 3.0 None None]
 [None None 3.0 4.0 5.0 5.0 4.0 3.0 None None]
 [None None 3.0 4.0 4.0 4.0 4.0 3.0 None None]
 [None None 3.0 3.0 3.0 3.0 3.0 3.0 None None]
 [None None None None None None None None None None]
 [None None None None None None None None None None]]

It seems that the searchlight did not run at the edge of the input data.
I try padding zero at the edge of data and run the searchlight again.

rad = 3 # specify searchlight radius
# creating a zero matrix with expanded x, y, z axis
paddata = np.zeros((dsize+rad*2, dsize+rad*2, dsize+rad*2, t))
# filling the original data at the center
paddata[rad:(dsize+rad), rad:(dsize+rad), rad:(dsize+rad), :] = data
padmask = np.zeros((dsize+rad*2, dsize+rad*2, dsize+rad*2))
padmask[rad:(dsize+rad), rad:(dsize+rad), rad:(dsize+rad)] = mask
print(f"padding 0 data shape: {paddata.shape}")
print(paddata[(dsize+rad*2)//2, :, :, 0])
print(f"padding 0 mask shape: {padmask.shape}")
print(padmask[(dsize+rad*2)//2, :, :])

padding 0 data shape: (16, 16, 16, 20)
[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0. 0. 0.]
 [0. 0. 0. 1. 2. 2. 2. 2. 2. 2. 2. 2. 1. 0. 0. 0.]
 [0. 0. 0. 1. 2. 3. 3. 3. 3. 3. 3. 2. 1. 0. 0. 0.]
 [0. 0. 0. 1. 2. 3. 4. 4. 4. 4. 3. 2. 1. 0. 0. 0.]
 [0. 0. 0. 1. 2. 3. 4. 5. 5. 4. 3. 2. 1. 0. 0. 0.]
 [0. 0. 0. 1. 2. 3. 4. 5. 5. 4. 3. 2. 1. 0. 0. 0.]
 [0. 0. 0. 1. 2. 3. 4. 4. 4. 4. 3. 2. 1. 0. 0. 0.]
 [0. 0. 0. 1. 2. 3. 3. 3. 3. 3. 3. 2. 1. 0. 0. 0.]
 [0. 0. 0. 1. 2. 2. 2. 2. 2. 2. 2. 2. 1. 0. 0. 0.]
 [0. 0. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]
padding 0 mask shape: (16, 16, 16)
[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0. 0. 0.]
 [0. 0. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0. 0. 0.]
 [0. 0. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0. 0. 0.]
 [0. 0. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0. 0. 0.]
 [0. 0. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0. 0. 0.]
 [0. 0. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0. 0. 0.]
 [0. 0. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0. 0. 0.]
 [0. 0. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0. 0. 0.]
 [0. 0. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0. 0. 0.]
 [0. 0. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]

# run the searchlight again with padding data
sl = Searchlight(sl_rad = rad, max_blk_edge=10)
sl.distribute([paddata], padmask)
sl.broadcast(0)
sl_results = sl.run_searchlight(sl_func, pool_size = 1)
print(f"output shape: {sl_results.shape}")
print(sl_results[(dsize+rad*2)//2, :, :])

output shape: (16, 16, 16)
[[None None None None None None None None None None None None None None None None]
 [None None None None None None None None None None None None None None None None]
 [None None None None None None None None None None None None None None None None]
 [None None None 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 None None None]
 [None None None 1.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 1.0 None None None]
 [None None None 1.0 2.0 3.0 3.0 3.0 3.0 3.0 3.0 2.0 1.0 None None None]
 [None None None 1.0 2.0 3.0 4.0 4.0 4.0 4.0 3.0 2.0 1.0 None None None]
 [None None None 1.0 2.0 3.0 4.0 5.0 5.0 4.0 3.0 2.0 1.0 None None None]
 [None None None 1.0 2.0 3.0 4.0 5.0 5.0 4.0 3.0 2.0 1.0 None None None]
 [None None None 1.0 2.0 3.0 4.0 4.0 4.0 4.0 3.0 2.0 1.0 None None None]
 [None None None 1.0 2.0 3.0 3.0 3.0 3.0 3.0 3.0 2.0 1.0 None None None]
 [None None None 1.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 2.0 1.0 None None None]
 [None None None 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 None None None]
 [None None None None None None None None None None None None None None None None]
 [None None None None None None None None None None None None None None None None]
 [None None None None None None None None None None None None None None None None]]

Then I could trim down the padding value and got more plausible results.

[mihaic]

Hello, @yushiangsu. Thank you for thoroughly documenting the issues. This is indeed the intended behavior: the outer sl_rad voxels are considered padding. We could change the implementation to add and remove the padding automatically, but that would be a backward-incompatible change. Instead, I think we should document the current behavior.