import numpy as np
from scipy.ndimage import median
 
dat = np.arange(12).reshape(2, 2, 3)
idx = np.array([[0, 0], [1, 2]])
 
def summarize(dat, idx):
    idx = np.unique(idx, return_inverse=True)[1].reshape(idx.shape)
    chan = dat.shape[-1]
    offset = idx.max() + 1
    index = np.stack([idx + i * offset for i in range(chan)], axis=-1)
    return median(dat, index, index=range(offset * chan)).reshape(chan, offset).T
 
print(summarize(dat, idx))

aW1wb3J0IG51bXB5IGFzIG5wCmZyb20gc2NpcHkubmRpbWFnZSBpbXBvcnQgbWVkaWFuCgpkYXQgPSBucC5hcmFuZ2UoMTIpLnJlc2hhcGUoMiwgMiwgMykKaWR4ID0gbnAuYXJyYXkoW1swLCAwXSwgWzEsIDJdXSkKCmRlZiBzdW1tYXJpemUoZGF0LCBpZHgpOgogICAgaWR4ID0gbnAudW5pcXVlKGlkeCwgcmV0dXJuX2ludmVyc2U9VHJ1ZSlbMV0ucmVzaGFwZShpZHguc2hhcGUpCiAgICBjaGFuID0gZGF0LnNoYXBlWy0xXQogICAgb2Zmc2V0ID0gaWR4Lm1heCgpICsgMQogICAgaW5kZXggPSBucC5zdGFjayhbaWR4ICsgaSAqIG9mZnNldCBmb3IgaSBpbiByYW5nZShjaGFuKV0sIGF4aXM9LTEpCiAgICByZXR1cm4gbWVkaWFuKGRhdCwgaW5kZXgsIGluZGV4PXJhbmdlKG9mZnNldCAqIGNoYW4pKS5yZXNoYXBlKGNoYW4sIG9mZnNldCkuVAoKcHJpbnQoc3VtbWFyaXplKGRhdCwgaWR4KSk=