I am trying to compute the histogram of an ND image. I have been able to use the array operation for a part of it and that makes it quite fast. There is one part of my code that I have not been able to vectorize and is therefore quite slow. I was wondering if there is a faster way to do this, since this is substantially slower than what I would get from C code.
# a and b are arrays (ND)
sz = product(a.shape) ai = a.flat bi = b.flat
for i in range(sz): hist_j[bi[i], ai[i]] += 1 hist_a[ai[i]] += 1 hist_b[bi[i]] += 1