On Wed, Jul 16, 2008 at 17:12, <Jack.Cook@shell.com> wrote:
Robert,
I can understand how this works if K is a constant time value but in my case K varies at each location in the two-dimensional slice. In other words, if I was doing this in a for loop I would do something like this
for i in range(numI): for j in range(numJ): k = slice(i,j) trace = cube(i,j,k-half_width:k+half_width) # shove trace in sub volume
What am I missing?
Ah, okay. It's a bit tricky, though. Yes, you need to use fancy indexing. Since axis you want to be index fancifully is not the first one, you have to be more explicit than you might otherwise want. For example, it would be great if you could just use slices for the first two axes: cube[:,:,slice + numpy.arange(-half_width,half_width)] but the semantics of that are a bit different for reasons I can explain later, if you want. Instead, you have to have explicit fancy-index arrays for the first two axes. Further, the arrays for each axis need to be broadcastable to each other. Fancy indexing will iterate over these broadcasted arrays in parallel with each other to form the new array. The liberal application of numpy.newaxis will help us achieve that. So this is the complete recipe: In [29]: import numpy In [30]: ni, nj, nk = (10, 15, 20) # Make a fake data cube such that cube[i,j,k] == k for all i,j,k. In [31]: cube = numpy.empty((ni,nj,nk), dtype=int) In [32]: cube[:,:,:] = numpy.arange(nk)[numpy.newaxis,numpy.newaxis,:] # Pick out a random fake horizon in k. In [34]: kslice = numpy.random.randint(5, 15, size=(ni, nj)) In [35]: kslice Out[35]: array([[13, 14, 9, 12, 12, 11, 8, 14, 11, 13, 13, 13, 8, 11, 8], [ 7, 12, 12, 6, 10, 12, 9, 11, 13, 9, 14, 11, 5, 12, 12], [ 7, 5, 10, 9, 6, 5, 5, 14, 5, 6, 7, 10, 6, 10, 11], [ 6, 9, 11, 14, 7, 11, 10, 6, 6, 9, 9, 11, 5, 5, 14], [12, 8, 11, 6, 10, 8, 5, 9, 8, 10, 7, 5, 9, 9, 14], [ 9, 8, 10, 9, 10, 12, 10, 10, 6, 10, 11, 6, 8, 7, 7], [11, 12, 7, 13, 5, 5, 8, 14, 5, 14, 9, 10, 12, 7, 14], [ 7, 7, 7, 12, 10, 6, 13, 13, 11, 13, 8, 11, 13, 14, 14], [ 6, 13, 13, 10, 10, 14, 10, 8, 9, 14, 13, 12, 9, 9, 5], [13, 14, 10, 8, 11, 11, 10, 6, 12, 11, 12, 12, 13, 11, 7]]) In [36]: half_width = 3 # These two replace the empty slices for the first two axes. In [37]: idx_i = numpy.arange(ni)[:,numpy.newaxis,numpy.newaxis] In [38]: idx_j = numpy.arange(nj)[numpy.newaxis,:,numpy.newaxis] # This is the substantive part that actually picks out our window. In [41]: idx_k = kslice[:,:,numpy.newaxis] + numpy.arange(-half_width,half_width+1) In [42]: smallcube = cube[idx_i,idx_j,idx_k] In [43]: smallcube.shape Out[43]: (10, 15, 7) # Now verify that our window is centered on kslice everywhere: In [47]: smallcube[:,:,3] Out[47]: array([[13, 14, 9, 12, 12, 11, 8, 14, 11, 13, 13, 13, 8, 11, 8], [ 7, 12, 12, 6, 10, 12, 9, 11, 13, 9, 14, 11, 5, 12, 12], [ 7, 5, 10, 9, 6, 5, 5, 14, 5, 6, 7, 10, 6, 10, 11], [ 6, 9, 11, 14, 7, 11, 10, 6, 6, 9, 9, 11, 5, 5, 14], [12, 8, 11, 6, 10, 8, 5, 9, 8, 10, 7, 5, 9, 9, 14], [ 9, 8, 10, 9, 10, 12, 10, 10, 6, 10, 11, 6, 8, 7, 7], [11, 12, 7, 13, 5, 5, 8, 14, 5, 14, 9, 10, 12, 7, 14], [ 7, 7, 7, 12, 10, 6, 13, 13, 11, 13, 8, 11, 13, 14, 14], [ 6, 13, 13, 10, 10, 14, 10, 8, 9, 14, 13, 12, 9, 9, 5], [13, 14, 10, 8, 11, 11, 10, 6, 12, 11, 12, 12, 13, 11, 7]]) In [50]: (smallcube[:,:,3] == kslice).all() Out[50]: True Clear as mud? I can go into more detail if you like particularly about how newaxis works. -- Robert Kern "I have come to believe that the whole world is an enigma, a harmless enigma that is made terrible by our own mad attempt to interpret it as though it had an underlying truth." -- Umberto Eco