@ptrblck thank you for your quick response, it works perfectly! However, for a single channel image I had to adapt it a little:
idx1, idx2, idx3 = indices.chunk(3, dim=3)
images[idx1, idx2, idx3].squeeze(3)