function merge,images,bboxes,bbox=bbox,replace=replace,add=add,brightest=brightest,dimmest=dimmest,average=average,nozero=nozero,first=first,last=last,mgcount=mgcount

  w=where(ptr_valid(images),nw)
  if nw eq 0 then return,0

  frames=n_elements(images); This includes the null images
  type=size((*(images[w[0]])),/type);
  if n_elements(first) eq 0 then first=0
  if n_elements(last) eq 0 then last=frames-1

  bbox=[min(bboxes[w,0]),min(bboxes[w,1]),max(bboxes[w,2]),max(bboxes[w,3])]
  mg=make_array(bbox[2]-bbox[0]+1,bbox[3]-bbox[1]+1,type=type,value=!values.f_nan);
  mgCount=make_array(bbox[2]-bbox[0]+1,bbox[3]-bbox[1]+1,type=2,value=0)

  for i=first,last do begin
    if ptr_valid(images[i]) then begin
      thisimage=(*(images[i]))
      xminprojmg=bboxes[i,0]-bbox[0];
      xmaxprojmg=bboxes[i,2]-bbox[0];
      yminprojmg=bboxes[i,1]-bbox[1];
      ymaxprojmg=bboxes[i,3]-bbox[1];
      submg=Mg[xminprojmg:xmaxprojmg,yminprojmg:ymaxprojmg]
      submgCount=MgCount[xminprojmg:xmaxprojmg,yminprojmg:ymaxprojmg]
      if keyword_set(replace) then begin
        if keyword_set(nozero) then begin
          w=where(logical_and(thisimage gt 0,finite(thisimage)),wc)
        end else begin
          w=where(finite(thisimage),wc)
        end
        if(wc gt 0) then begin
          submgCount[w]+=1
          subMg[w]=thisimage[w]
        end
      end
      if keyword_set(brightest) or keyword_Set(dimmest) then begin
        if keyword_set(nozero) then begin
          wf=where(logical_and(finite(subMg),logical_and(thisimage gt 0,finite(thisimage))),wfc)
          wi=where(logical_and(~finite(subMg),logical_and(thisimage gt 0,finite(thisimage))),wic)

        end else begin
          wf=where(logical_and(finite(subMg),finite(thisimage)),wfc)
          wi=where(logical_and(~finite(subMg),finite(thisImage)),wic)
        end
        if(wic gt 0) then begin
          submgCount[wi]+=1
          submg[wi]=thisImage[wi]
        end
        if(wfc gt 0) then begin
          if keyword_set(brightest) then begin
            subMg[wf]=thisimage[wf]>subMg[wf]
          end else begin
            subMg[wf]=thisimage[wf]<subMg[wf]
          end
          subMgCount[wf]+=1
        end
      end
      if keyword_set(add) or keyword_set(average) then begin
        subMg=total([[[submg]],[[thisImage]]],3,/nan)
        if keyword_set(nozero) then begin
          w=where(logical_and(thisimage gt 0,finite(thisimage)),wc)
        end else begin
          w=where(finite(thisimage),wc)
        end
        if wc gt 0 then submgCount[w]+=1
      end

      Mg[xminprojmg:xmaxprojmg,yminprojmg:ymaxprojmg]=submg
      MgCount[xminprojmg:xmaxprojmg,yminprojmg:ymaxprojmg]=submgCount
    end
  endfor
  if keyword_set(average) then begin
    w=where(mgCount gt 0)
    mg[w]/=mgCount[w]
  end
  return,Mg
end