pro associate, l_assoc, labels
    ll = labels[uniq(labels, sort(labels))]
    n = N_ELEMENTS(ll)
    for i=0L,n-1 do $
        ll = [ll, l_assoc[ll[i]], where(l_assoc eq ll[i])]
    ll = ll[uniq(ll, sort(ll))]
    if ll[0] eq 0 then ll = ll[1:N_ELEMENTS(ll)-1]
    i = l_assoc[ll[0]] ; choose the least element (already sorted)
    ii = ll
    l_assoc[ii] = i
    ;print, 'labels: ', ll
    ;print, l_assoc
end

; calls label_region with gutters around the data
function safe_label_region, data
    s = size(data)
    d = s[1:s[0]]
    i = INDGEN(8) lt s[0]
    i2 = INTARR(8) & i2[0] = d - 1
    i2 += i
    data2 = BYTARR(d+2)
    data2[i[0],i[1],i[2],i[3],i[4],i[5],i[6],i[7]] = data ne 0
    return, (label_region(data2,/ALL_NEIGHBORS,/ULONG))[ $
        i[0]:i2[0],i[1]:i2[1],i[2]:i2[2],i[3]:i2[3],i[4]:i2[4], $
        i[5]:i2[5],i[6]:i2[6],i[7]:i2[7]]
end

pro ghost_labels_and_associate, label, dims, l_assoc
    ; order of label and dims is
    ; self parent [neighbors]
    
    s = dims[3:5,0]+1-dims[0:2,0]
    ; first, make the working array
    labels = ULONARR([s+2,2])
    ; splat the parent
    p_start = (( dims[0:2,0]    / 2)     - 1) - dims[0:2,1]
    p_end   = (((dims[3:5,0]+1) / 2) - 1 + 1) - dims[0:2,1]
    p_i0 = p_start > 0 < (dims[3:5,1]-dims[0:2,1])
    p_i1 = p_end   > 0 < (dims[3:5,1]-dims[0:2,1])
    
    p = ULONARR(s/2 + 2)
    p[p_i0[0]-p_start[0],p_i0[1]-p_start[1],p_i0[2]-p_start[2]] = $
        (*label[1])[p_i0[0]:p_i1[0],p_i0[1]:p_i1[1],p_i0[2]:p_i1[2]]
    p = REBIN(p, s+4, /SAMPLE)
    labels[0,0,0,0] = p[1:s[0]+2,1:s[1]+2,1:s[2]+2]
    labels[0,0,0,1] = labels[*,*,*,0]

    ; splat the neighbors
    if N_ELEMENTS(label) gt 2 then begin
        n = 2 + INDGEN(N_ELEMENTS(label)-2)
        left    = CMREPLICATE(dims[0:2,0]-1, N_ELEMENTS(n))
        right   = CMREPLICATE(dims[3:5,0]+1, N_ELEMENTS(n))
        n_left  = dims[0:2,n] > left < right
        n_right = dims[0:2,n] > left < right
        ; nn_s: neighbor's neighbor-index coords, start
        ; ns_e: neighbor's splat-index coords, end
        nn_s = n_left  - dims[0:2,n]
        nn_e = n_right - dims[0:2,n]
        ns_s = n_left  - left
        ns_e = n_right - left
        for i=0L,N_ELEMENTS(n)-1 do begin
            labels[ns_s[0,i]:ns_e[0,i], ns_s[1,i]:ns_e[1,i], ns_s[2,i]:ns_e[2,i], 0] = $
                (*label[n[i]])[nn_s[0,i]:nn_e[0,i], nn_s[1,i]:nn_e[1,i], nn_s[2,i]:nn_e[2,i]]
        endfor
    endif
    
    ; splat the grid
    labels[1,1,1,0] = *label[0]

    ; now generate an association map
    l = safe_label_region(labels gt 0)
    h = histogram(l, REVERSE_INDICES=ri, BINSIZE=1, MIN=0)

    for i=1L,N_ELEMENTS(h)-1 do begin
        ; associate all of the labels in this region
        ll = labels[ri[ri[i]:ri[i+1]-1]]
        ll = l_assoc[ll[uniq(ll,sort(ll))]]
        if N_ELEMENTS(ll) gt 1 then $
            associate, l_assoc, ll
    endfor
end

; segmentation data structure:
; {segmentation, n_labels:LONG, labels_index:PTRARR(LONARR)}
; array structure:
; [0:n_grids-1] self-index of first grid-index of this label for grid i
; [n_grids]     self-index pointing after last grid-index

function segment_volume, PROJ=proj, threshold, HISTOGRAM=plot_histogram, BINSIZE=bin
    ; Fetch data
        if N_ELEMENTS(proj) eq 0 then $
            proj = build_projection(INTERPOLATION=0, /NOSWAPPING)
        grids = *proj.grids
        dims = grid_int_coords(INFO=grids)
        neighbors = grid_neighbors(INFO=grids, DIMS=dims, CELLS=1)
        get_projection_data, PROJECTION=proj, FIELDS=["Density"], /NOPARTS
        grid_data = *proj.grid_data
    
    ; Threshold and classify
        label = PTRARR(N_ELEMENTS(grid_data), /ALLOCATE_HEAP)
        next_region = 0L
        for j=0L,N_ELEMENTS(grid_data)-1 do begin
            m = *grid_data[j].data[0] * proj.density_unit
            for k=0L,N_ELEMENTS(grid_data)-1 do begin
                if grids[k].parent eq grids[j].num then begin
                    i0 = grids[k].parent_s_index
                    i1 = grids[k].parent_e_index
                    m[i0[0]:i1[0],i0[1]:i1[1],i0[2]:i1[2]] = 0.d
                endif
            endfor
            l = safe_label_region(m ge (threshold * 2.33)) ; mean molecular weight
            ; uniquely number regions
            n_regions = MAX(l)
            ind = where(l ne 0)
            if ind[0] ge 0 then $
                l[ind] += next_region
            *label[j] = l
            next_region += n_regions
        endfor
    
    ; Find associations
        l_assoc = ULINDGEN(next_region+1) ; 1-based
        order = REVERSE(SORT(grids.level))
        for i=0L,N_ELEMENTS(grid_data)-1 do begin
            j = order[i]
            p = (where(grids.num eq grids[j].parent))[0]
            if p lt 0 then continue
            ind = [j,p]
            n_ind = where(neighbors[*,j] ge 0)
            if n_ind[0] ge 0 then $
                ind = [ind,neighbors[n_ind, j]]
            ghost_labels_and_associate, label[ind], dims[*,ind], l_assoc
        endfor
        ;print, l_assoc

    ; Build output index lists
        labels = l_assoc[1:N_ELEMENTS(l_assoc)-1]
        orig_label = labels[uniq(labels, sort(labels))]
        for i=0L,N_ELEMENTS(orig_label)-1 do l_assoc[where(l_assoc eq orig_label[i])] = i+1
        n_labels = N_ELEMENTS(orig_label)
        n_grids = N_ELEMENTS(grid_data)
        segmentation = PTRARR(n_labels, /ALLOCATE_HEAP)
        index_count = LON64ARR(n_labels, n_grids)
        for j=0L,n_grids-1 do begin
            l = l_assoc[*label[j]]
            index_count[0,j] = histogram(l, BINSIZE=1, MIN=1, MAX=n_labels)
            *label[j] = l
        endfor
        for i=0L,n_labels-1 do begin
            *segmentation[i] = [n_grids+1,TOTAL(REFORM(index_count[i,*]),/CUMULATIVE,/INTEGER)+n_grids+1,LONARR(TOTAL(index_count[i,*]))]
        endfor
        for j=0L,N_ELEMENTS(grid_data)-1 do begin
            l = *label[j]
            h = histogram(l, REVERSE_INDICES=r, BINSIZE=1, MIN=1, MAX=n_labels)
            for i=0L,n_labels-1 do begin
                if r[i+1] eq r[i] then continue
                (*segmentation[i])[(*segmentation[i])[j]] = r[r[i]:r[i+1]-1]
            endfor
        endfor
    return, {segmentation,n_labels:n_labels,labels_index:PTR_NEW(segmentation)}
end

pro segmentation_statistics, PROJ=proj, SEGMENTS=seg
    if N_ELEMENTS(proj) eq 0 then $
        proj = build_projection(INTERPOLATION=0, /NOSWAPPING)
    grids = *proj.grids
    dims = grid_int_coords(INFO=grids)
    neighbors = grid_neighbors(INFO=grids, DIMS=dims, CELLS=1)
    get_projection_data, PROJECTION=proj, FIELDS=["Density", "x-velocity", "y-velocity", "z-velocity"], /NOPARTS
    grid_data = *proj.grid_data
    rho_i = (where(grid_data[0].data_fields eq "Density"))[0]
    vx_i  = (where(grid_data[0].data_fields eq "x-velocity"))[0]
    vy_i  = (where(grid_data[0].data_fields eq "y-velocity"))[0]
    vz_i  = (where(grid_data[0].data_fields eq "z-velocity"))[0]

    protonmass = 1.67e-24
    dV = REFORM(PRODUCT(*proj.delta_distance, 1)) * proj.distance_unit^3
    dp = (*proj.delta_distance)[2,*]

    n_cores = seg.n_labels
    n_grids = N_ELEMENTS(grids)
    core_mass = DBLARR(n_cores)
    core_volume = DBLARR(n_cores)
    core_KE = DBLARR(n_cores)
    core_pos = DBLARR(n_cores, 3)
    core_bounds = DBLARR(n_cores, 3, 2)
    core_bounds[*,*,0] = 1.0

    image = DBLARR(proj.slice.image_size, proj.slice.image_size)
    labels_index = *seg.labels_index
    ri = LON64ARR(n_grids+1, n_cores)
    for i=0L,n_cores-1 do ri[0,i] = (*labels_index[i])[0:n_grids]

    for j=0L,N_ELEMENTS(grid_data)-1 do begin
        m  = *grid_data[j].data[rho_i] * proj.density_unit * protonmass / 1.98892d33 ; solar masses
        vx = *grid_data[j].data[vx_i]
        vy = *grid_data[j].data[vy_i]
        vz = *grid_data[j].data[vz_i]
        d = TRANSPOSE(CMREPLICATE((size(m))[1:3],N_ELEMENTS(m))) ; cells per column (1), cells per row (x), cells per plane (x*y)
        s = TRANSPOSE(CMREPLICATE([1, PRODUCT((size(m))[1:2], /CUMULATIVE, /INTEGER)],N_ELEMENTS(m))) ; cells per column (1), cells per row (x), cells per plane (x*y)
        xyz = (CMREPLICATE(TOTAL(dims[0:2,j] * s, /INTEGER) + LONARR(N_ELEMENTS(m)), 3) / s) mod d
        xyz = ISHFT(TEMPORARY(xyz), max(grids.level) - grids[j].level)
        for i=0L,n_cores-1 do begin
            ki = ri[j:j+1, i]
            if ki[0] eq ki[1] then continue
            k = (*labels_index[i])[ki[0]:ki[1]-1]

            core_mass[i]      += TOTAL(m[k]) * dV[j]
            core_volume[i]    += N_ELEMENTS(k) * dV[j]
            core_KE[i]        += .5 * TOTAL(m[k] * (vx[k]^2+vy[k]^2+vz[k]^2)) * dV[j]
            core_pos[i,0]      = core_pos[i,*] + TOTAL(CMREPLICATE(m[k] * dV[j], 3) * xyz[k,*], 1)
            core_bounds[i,0,0] = core_bounds[i,*,0] < min(xyz[k,*],DIMENSION=1)
            core_bounds[i,0,1] = core_bounds[i,*,1] > max(xyz[k,*],DIMENSION=1)
        endfor
    endfor
    core_pos /= CMREPLICATE(core_mass, 3)
    for i=0L,n_cores-1 do begin
        print, "Core ", i+1, ": m=", core_mass[i], ", l=", core_volume[i]^.3, ", pos=[",core_pos[i,0],",",core_pos[i,1],",",core_pos[i,2],"]"
    endfor
    bin = 0.2
    h = histogram(alog10(core_mass), BINSIZE=bin, OMIN=omin, OMAX=omax)
    x = 10^(OMIN+DINDGEN(N_ELEMENTS(h))*bin)
    plot, x, h, PSYM=10, /XLOG, /YLOG, YRANGE=[1,max(h)]
    m = max(h,i)
    oplot, x, m*(x^(-1.35))/(x[i]^(-1.35))
end

;pro render_segmentation
;    ;grid_xy = *proj.grid_xy
;    ;grid_wh = *proj.grid_wh
;            ;xy = grid_xy[0:1,j]
;            ;wh = grid_wh[*,j]
;            ;dst = [xy,xy+wh-1] > 0 < (proj.slice.image_size-1)
;            ;src = dst - [xy,xy]
;            ;image[dst[0]:dst[2], dst[1]:dst[3]] += (CONGRID(TOTAL(m,3)*dp[j], dst[2]-dst[0]+1, dst[3]-dst[1]+1, /CENTER))[src[0]:src[2],src[1]:src[3],*]
;    if not KEYWORD_SET(plot_histogram) then begin 
;        tvscl, (image)
;        ; exclude cores with less than 1 msol
;        ind = where(core_mass gt 0.5d)
;        if ind[0] ge 0 then begin
;            w = ISHFT(dims[3:4,0] - dims[0:1,0] + 1, max(grids.level)-grids[0].level)
;            plot, core_pos[ind,0], core_pos[ind,1], PSYM=5,XMARGIN=[0,0],YMARGIN=[0,0],/NOERASE,XRANGE=[0,w[0]],YRANGE=[0,w[1]]
;        endif
;    endif else begin
;        h = histogram(alog10(core_mass), BINSIZE=bin, OMIN=omin, OMAX=omax)
;        x = 10^(OMIN+DINDGEN(N_ELEMENTS(h))*bin)
;        plot, x, h, PSYM=10, /XLOG, /YLOG, YRANGE=[1,max(h)]
;        m = max(h,i)
;        oplot, x, m*(x^(-1.35))/(x[i]^(-1.35))
;    endelse
;end