function define_slice, SLICE_ORIENTATION=ori, SLICE_SIZE=xy_size, SLICE_VALUE=z_center, DEPTH_VALUE=z_range, IMAGE_SIZE=image_size
    ori = DOUBLE(ori)
    Z_axis = ori / NORM(ori)

    ; which axis should we project along?
    z_rate = max(abs(Z_axis), Z) ; rate at which 'photons' progress along the projection axis
    Z += 3*(Z_axis[Z] lt 0) ; determine whether Z is reversed
    ; which zone does this correspond to?
    n = ([0,1,4,2,3,5])[Z]
    ; determine the spherical angles
    alpha = atan(ori[1], ori[0])
    beta  = atan(sqrt(ori[0]^2+ori[1]^2), ori[2])
    ; now look up the permutation and coefficients
    ca = cos(alpha) & sa = sin(alpha) & ta = tan(alpha)
    cb = cos(beta ) & sb = sin(beta ) & tb = tan(beta )
    coeffs = ([[4,2,0, ca,  0, sa*cb,    sb,    ta           ,-(sa*ta+ca)/tb        ], $
               [0,2,1, sa,  0,-ca*cb,    sb,   -1d/ta        ,-(ca/ta+sa)/tb        ], $
               [1,2,3,-ca,  0,-sa*cb,    sb,    ta           , (sa*ta+ca)/tb        ], $
               [3,2,4,-sa,  0, ca*cb,    sb,   -1d/ta        , (ca/ta+sa)/tb        ], $
               [1,0,2,-ca, sa,-sa*cb,-ca*cb,-ta*tb/(sa*ta+ca),        -tb/(sa*ta+ca)], $
               [1,3,5,-ca,-sa,-sa*cb, ca*cb, ta*tb/(sa*ta+ca),        -tb/(sa*ta+ca)]])[*,n]
    x = FIX(coeffs[0]) & y = FIX(coeffs[1]) & z = FIX(coeffs[2])
    a = coeffs[3] & b = coeffs[4] & c = coeffs[5] & d = coeffs[6] & su = coeffs[7] & sv = coeffs[8]

    xi = x mod 3 & yi = y mod 3 & zi = z mod 3 & xyi = [xi,yi]
    axes = [X,Y,Z]

    bounds_data = DBLARR(3,2)
    bounds_proj = DBLARR(3,2)

    ; bounds in the projection plane
    ; adjust xy_size ever so slightly to make sure the left edge coincides with
    ; an integer multiple of the finest resolution pixel width
    dxdy = DOUBLE(xy_size[2:3]-xy_size[0:1])/image_size
    bounds_proj[0:1,0] = dxdy*LONG64(DOUBLE(xy_size[0:1])/dxdy[0:1])
    bounds_proj[0:1,1] = bounds_proj[0:1,0]+dxdy*image_size

    ; bounds along the axis of projection
    z_max = z_center+z_range/2.  ;z_max = nearest_valid_depth_value(z_max, grid_info[0].dim[0]) < 1
    z_min = z_max-z_range
    if z_min lt 0. then z_max = z_range
    if z_max gt 1. then z_min = 1.-z_range
    bounds_proj[2,*] = [[z_min], [z_max]]
    
    bounds_data[[xi,yi,zi,xi+3,yi+3,zi+3]] = bounds_proj
    if axes[0] ge 3 then bounds_data[xi,*] = 1.d - bounds_data[xi,[1,0]]
    if axes[1] ge 3 then bounds_data[yi,*] = 1.d - bounds_data[yi,[1,0]]
    if axes[2] ge 3 then bounds_data[zi,*] = 1.d - bounds_data[zi,[1,0]]
    return, {slice, axes:axes, bounds_proj:bounds_proj, bounds_data:bounds_data, xyi:xyi, xi:xi, yi:yi, zi:zi, image_size:image_size, $
             alpha:alpha, beta:beta, z_rate:z_rate, a:a, b:b, c:c, d:d, su:su, sv:sv}
end