function composite_HII_picture, im, LEVELS=levels, DRY_RUN=dry_run
    im1 = im[*,*,0] & im2 = im[*,*,1] & im3 = im[*,*,2] & im4 = im[*,*,3]

    im3 = im3^.3
    blob_profile = [1.,.3,.1,.03,.01,.003]
    blob = blob_profile[[[5,4,3,4,5],[4,2,1,2,4],[3,1,0,1,3],[4,2,1,2,4],[5,4,3,4,5]]]
    im3 = convol(im3, blob, MISSING=0., /EDGE_ZERO)
    
    read_gif, 'star.gif', star_data, rr, gg, bb
    star_data = (star_data/255.0)^3
    star_data /= total(star_data) ; normalize
    im4 = convol(im4^.01, star_data, MISSING=0., /EDGE_ZERO)

    if KEYWORD_SET(dry_run) or (N_ELEMENTS(levels) eq 0) then begin
        levels = DBLARR(4)

        smoothed = median(im1, 5)
        levels[0] = (mean(smoothed) + 7*stddev(smoothed))
        smoothed = median(im2, 5)
        levels[1] = (mean(smoothed) + 7*stddev(smoothed))

        ind = where(im3 gt 0.)
        if ind[0] ge 0 then levels[2] = (mean(im3[ind]) + 2*stddev(im3[ind])) $
        else levels[2] = 1.D

        ind = where(im4 gt 0.)
        if ind[0] ge 0 then levels[3] = mean(im4[ind])*0.001 $
        else levels[3] = 1.D
    endif
    if not KEYWORD_SET(dry_run) then begin
        im1 = (im1 < levels[0]) / levels[0]
        im2 = (im2 < levels[1]) / levels[1]
        im3 = (im3 < levels[2]) / levels[2]
        im4 = (im4 < levels[3]) / levels[3]

        im5 = 1.D - (1.D - im3)*(1.D - im4)

        im = transpose([[[im1^.3]],[[im2^.6]],[[im2^.6]]], [2,0,1])
        im[1,*,*] = .8*im[1,*,*]+.2*(im[0,*,*])
        im = im > 0 < 1
        im[0,*,*] = 1 - (1-im[0,*,*]) * (1 - im5)^2
        im[1,*,*] = 1 - (1-im[1,*,*]) * (1 - im5)^2
        im[2,*,*] = 1 - (1-im[2,*,*]) * (1 - im5)^2
        return, BYTE(255. * im)
    endif else return, 1
end

function make_HII_picture, PROJ=proj, NO_NORMALIZE=no_normalize, VERBOSE=verbose, CACHE=cache
    if not KEYWORD_SET(proj) then $
        proj = build_projection(INTERP=1)
    no_normalize = KEYWORD_SET(no_normalize)
    if not KEYWORD_SET(cache) then begin
        cache = PTRARR(3, /ALLOCATE_HEAP)
    endif

    line_index = (where((line_emissivity(/AVAILABLE_LINES)).name eq 'H-alpha'))[0]
    if KEYWORD_SET(verbose) then set_spectral_palette, (line_emissivity(/AVAILABLE_LINES))[line_index].lambda
    im1 = project_with_absorption(PROJECTION=proj, LINE_INDEX=line_index, OPACITY=0, STARS_ONLY=0, VERB=verbose, CACHE=*cache[0])
    
    line_index = (where((line_emissivity(/AVAILABLE_LINES)).name eq '[OIII] 500.7nm'))[0]
    if KEYWORD_SET(verbose) then set_spectral_palette, (line_emissivity(/AVAILABLE_LINES))[line_index].lambda
    im2 = project_with_absorption(PROJECTION=proj, LINE_INDEX=line_index, OPACITY=0, STARS_ONLY=0, VERB=verbose, CACHE=*cache[1])

    if KEYWORD_SET(verbose) then loadct, 0
    im3 = project_with_absorption(PROJECTION=proj, LINE_INDEX=line_index, OPACITY=0, STARS_ONLY=1, VERB=verbose, CACHE=*cache[2])
    im4 = project_with_absorption(PROJECTION=proj, LINE_INDEX=line_index, OPACITY=0, STARS_ONLY=2, VERB=verbose, CACHE=*cache[2])
    im = [[[im1]],[[im2]],[[im3]],[[im4]]]
    if KEYWORD_SET(no_normalize) then $
        return, im $
    else $
        return, composite_HII_picture(im)
end