I've just written this snippet of code: result = numpy.empty((dim, dim, dim), numpy.bool) for x in xrange(dim): for y in xrange(dim): for z in xrange(dim): result[x, y, z] = ((xnear[y, z] < xfar[y, z]) and (ynear[x, z] < yfar[x, z]) and (znear[x, y] < zfar[x, y])) Is there a way to get numpy to do the looping?