function define_projection, GRIDS=grids, SLICE=slice, FASTMODE=min_pixel
    if not KEYWORD_SET(layers) then layers = 1
    if N_ELEMENTS(grids) lt 1 then begin
        print, 'define_projection: no grid information in grids:', grids
        STOP
    endif

    xi = slice.xi & yi = slice.yi & zi = slice.zi & xyi = slice.xyi
    ; use only grids that will show up on the image

    grid_flag       = grid_in_bounds(grids, slice.bounds_data)
    index_range     = grids.end_index[xi]-grids.start_index[xi]+1
    delta_distance  = DOUBLE(grids.right_edge[xi] - grids.left_edge[xi])/DOUBLE(index_range)
    del_pix         = DOUBLE(slice.image_size)*delta_distance/(slice.bounds_data[xi,1]-slice.bounds_data[xi,0])
    del_pix[0] = 1 & grid_flag[0] = 1 ; we always want at least the coarsest grid
    used_grids      = grids[where((del_pix ge min_pixel) AND (grid_flag gt 0))]
    n_grids         = N_ELEMENTS(used_grids)

    index_range     = used_grids.end_index-used_grids.start_index+1
    delta_distance  = (used_grids.right_edge - used_grids.left_edge)/DOUBLE(index_range)
    n_grids         = N_ELEMENTS(used_grids)
    ii_l            = TRANSPOSE((FLOOR((CMREPLICATE(slice.bounds_data[*,0], n_grids) - used_grids.left_edge) / delta_distance)    ) > 0 < (index_range-1))
    ii_r            = TRANSPOSE(( CEIL((CMREPLICATE(slice.bounds_data[*,1], n_grids) - used_grids.left_edge) / delta_distance) - 1) > 0 < (index_range-1))
    ii_d            = ii_r - ii_l + 1
    dp              = DOUBLE(slice.image_size)/CMREPLICATE(slice.bounds_data[xyi,1]-slice.bounds_data[xyi,0], N_ELEMENTS(used_grids)) * delta_distance

    ; get the dimensions of all of the grids in absolute index coords
    dims = grid_int_coords(INFO=used_grids)
    ; rescale to a common resolution and bound by ii_l, ii_r
    scale = TRANSPOSE(CMREPLICATE(max(used_grids.level) - [used_grids.level], 6))
    dims = ISHFT([dims[0:2,*]+TRANSPOSE(ii_l), dims[0:2,*]+TRANSPOSE(ii_r)+1], scale)
    ; find the number of maximum-resolution indices (equals # of coarse cells * 2^max level)
    int_width = ISHFT((used_grids[(where([used_grids.level] eq 0))[0]].dim), max(used_grids.level))
    ; convert to projection space by permuting axes
    dims[0,0] = dims[[xi,yi,zi,xi+3,yi+3,zi+3],*]
    iw = int_width
    if slice.axes[0] ge 3 then dims[[0,3],*] = iw[0] - dims[[3,0],*]
    if slice.axes[1] ge 3 then dims[[1,4],*] = iw[1] - dims[[4,1],*]
    if slice.axes[2] ge 3 then dims[[2,5],*] = iw[2] - dims[[5,2],*]
    ; find the dimensions of the slice in these units
    int_slice = REFORM(slice.bounds_proj, 6) * [iw,iw]

    ; ########## COMPUTE DEPTH LIST ########################################
    ; get cell grid points and sort them by depth. Remember which grids
    ; they correspond to at the same time. 
    count_total = LONG([0L, TOTAL(ii_d[*,zi], /CUMULATIVE)]) ; cumulative sum with 0 tacked on the front
    n_depths = count_total[N_ELEMENTS(used_grids)]
    int_depth_values = LON64ARR(n_depths) ; pre-allocate arrays
    depth_values = DBLARR(n_depths)
    grid_index = INTARR(n_depths)
    grid_values = LONARR(n_depths)
    for i=0L, n_elements(used_grids)-1 do begin
        int_depth_values[count_total[i]] = dims[2,i] + ISHFT(L64INDGEN(ii_d[i,zi]), scale[zi,i])
        grid_index[count_total[i]] = INDGEN(ii_d[i,zi])
        grid_values[count_total[i]] = LONARR(ii_d[i,zi])+i
    endfor
    unique_depths = UNIQ(int_depth_values, REVERSE(SORT(int_depth_values))) ; calculate index list for back-to-front iteration of unique depths
    depth_values = DOUBLE(int_depth_values) / int_width[zi]
    n_depths = N_ELEMENTS(unique_depths)

    int_slice_xy = int_slice[0:1]
    int_slice_wh = int_slice[3:4] - int_slice[0:1]
    ; compute x,y coordinates on the screen
    grid_xy = LONG((((dims[[0,1,3,4],*]-CMREPLICATE([int_slice_xy, int_slice_xy], n_grids)) * slice.image_size) / CMREPLICATE([int_slice_wh,int_slice_wh], n_grids)))
    grid_wh = grid_xy[2:3,*] - grid_xy[0:1,*] + 1 ; is this plus 1 correct? not sure

    distance_unit  = get_unit('Length', /FORCE_USE_UNITS) * 3.08568e18 ; centimeters
    intensity_unit = get_unit('Intensity', /FORCE_USE_UNITS)
    density_unit   = get_unit('Density', /FORCE_USE_UNITS)
    
    return, {projection, slice:slice, fieldnames:STRARR(100), $
            grids:PTR_NEW(used_grids), ii_l:PTR_NEW(ii_l), ii_r:PTR_NEW(ii_r), ii_d:PTR_NEW(ii_d), delta_distance:PTR_NEW(delta_distance), dp:PTR_NEW(dp), $
            depth_values:PTR_NEW(depth_values), grid_index:PTR_NEW(grid_index), grid_values:PTR_NEW(grid_values), unique_depths:PTR_NEW(unique_depths), n_depths:n_depths, $
            distance_unit:distance_unit, intensity_unit:intensity_unit, border:0L, density_unit:density_unit, $
            grid_xy:PTR_NEW(grid_xy), grid_wh:PTR_NEW(grid_wh), dims:PTR_NEW(dims), iw:iw, $
            grid_data:PTR_NEW(), part_data:PTR_NEW()}
end