function average_images, raw_mosaic, total_weight, imxmin, imymin, $
				     BORDER=border, BACK_SET=back, $
				     GROUP=group, WHICH_IMAGES=wim, $
				     FRAC_PIX=frac, WEIGHT_IMAGES=weights
;+
; NAME:
;	average_images
; PURPOSE:
;	Average the images in raw_mosaic and return the averaged_mosaic image.
; CALLING:
;	mosaic = average_images( raw_mosaic )
; INPUTS:
;	raw_mosaic = structured array of images with relative Location info.
; KEYWORDS:
;      /BACK_SET : for pixels with no data, set the value < minimum of data,
;		This allows them to be distinguished from all other data pixels.
;			Default is to leave them = 0.
;      /FRAC_PIX : if relative offsets between images are in fractional pixels,
;			use this info in creating final average mosaic.
;			Default is to round off to nearest pixel offset.
;	BORDER = width of border of each image to ignore (exclude from average),
;			default = 0.
;	GROUP = average only image in specified group number (default is all).
;	WHICH_IMAGES = the indices of images in raw_mosaic to average (def=all).
;	WEIGHT_IMAGES = weighting to use for each image when averaging,
;			default = 1 for all image.
; OUTPUTS:
;	total_weight = an image of the sum of weight in each pixel of
;			the final averaged mosaic. Thus if weights = 1,
;			each element of total_weight is the count of
;			how many images overlapped at that pixel.
;	imxmin, imymin = the location of bottom left corner of final image,
;			in relative offset coordinate.
; RESULT:
;	Function returns image with each pixel = average of overlapping pixels,
;	and this image is converted to the same type as the raw images.
; EXTERNAL CALLS:
;	function frac_pix_shift
;	function round_off
;	function conv_vartype
; PROCEDURE:
;	Loop thru images in raw_mosaic and add into big mosaic array,
;	also accumulating the total_weight array, then divide by it.
; HISTORY:
;	Written, Frank Varosi NASA/GSFC 1989.
;	F.V.1989, added keywords for BORDER = pixel border to ignore,
;				GROUP=group # to average,
;				/BACK_SET to set background empty pixels.
;	F.V.1991, optimized and added conv_vartype at end.
;	F.V.1991, added keyword WHICH_IMS= subscripts of images.
;	F.V.1992, option /FRAC_PIX to use fractional pixel offsets
;				in averaging (func frac_pix_shift()).
;-
	Nim = N_struct( raw_mosaic )
	if (Nim LE 0) then return,bytarr(2,2)

	if N_elements( border ) NE 1 then border=0
	border = fix( border > 0 )

	sim = size( raw_mosaic(0).image )
	xsiz = sim(1)
	xsizb = xsiz - border
	xsizb1 = xsizb - 1
	xsizbb = xsizb - border
	xsizbb1 = xsizbb - 1

	ysiz = sim(2)
	ysizb = ysiz - border
	ysizb1 = ysizb - 1
	ysizbb = ysizb - border
	ysizbb1 = ysizbb - 1

	if N_elements( wim ) LE 0 then wim = indgen( Nim )
	Nim = N_elements( wim )

	if N_elements( group ) EQ 1 then begin
		if (group GE 0) then begin
			wim = where( [raw_mosaic.group] EQ group, Nim )
			if (Nim LE 0) then begin
				print," group #",strtrim( group, 2 )," is empty"
				return,bytarr(2,2)
			   endif
			print," using images in group #", strtrim( group, 2 )
		   endif
	   endif

	print,fix( Nim )," images to average in mosaic"
	if (border GT 0) then print," ignoring border of width=",byte( border )

	if (Nim LE 0) then return, bytarr(2,2)
	if (Nim LE 1) then $
		return, raw_mosaic(wim(0)).image( border:xsizb1, border:ysizb1 )

	xmin = raw_mosaic(wim).Locx
	xmax = xmin + xsiz
	ymin = raw_mosaic(wim).Locy
	ymax = ymin + ysiz

	xminmin = min( xmin )
	xmaxmax = max( xmax )
	yminmin = min( ymin )
	ymaxmax = max( ymax )

	imxsiz = round_off( xmaxmax - xminmin ) - 2*border
	imysiz = round_off( ymaxmax - yminmin ) - 2*border
	imxmin = xminmin + border
	imymin = yminmin + border
	mosaic_image = make_array( imxsiz, imysiz, /FLOAT )

	if N_elements( weights ) LT Nim then weights = replicate( 1, Nim )
	s = size( weights )
	total_weight = make_array( imxsiz, imysiz, TYPE=s(s(0)+1) )

	xoff = xmin - xminmin
	yoff = ymin - yminmin
	xoffr = round_off( xoff )
	yoffr = round_off( yoff )
	xLastr = xoffr + xsizbb1
	yLastr = yoffr + ysizbb1

	if keyword_set( frac ) then begin
		xshift = xoff - xoffr
		yshift = yoff - yoffr
	  endif else begin
		xshift = fltarr( Nim )
		yshift = fltarr( Nim )
	   endelse

	for i=0,Nim-1 do begin

		image = frac_pix_shift( raw_mosaic(wim(i)).image, xshift(i), $
							/RENORMAL, yshift(i) )
		if (border GT 0) then $
				image = image( border:xsizb1, border:ysizb1 )

		if (weights(i) NE 1) then image = weights(i) * image

		xo = xoffr(i)
		yo = yoffr(i) 
		xL = xLastr(i)
		yL = yLastr(i)

		mosaic_image(xo,yo) = image + mosaic_image( xo:xL, yo:yL )
				   
		total_weight(xo,yo) = total_weight( xo:xL, yo:yL ) + $
					replicate( weights(i), xsizbb, ysizbb )
		print,form="($,i3)",i
	  endfor

	npixtot = Long( imxsiz ) * imysiz
	print," "
	print, strtrim( imxsiz,2 )," x ", strtrim( imysiz,2 ), $
		" = ",strtrim( npixtot,2 )," total pixels in mosaic..."

	wz = where( total_weight LE 0, nzero )
	pczero = (100.*nzero)/npixtot
	print,nzero," empty pixels (", pczero, " % )"

	type_out = sim(sim(0)+1)
	if !DEBUG then stop

	if (nzero LE 0) then begin

		print," ...averaging..."
		return, conv_vartype( mosaic_image/total_weight, TYPE=type_out )

	 endif else if (pczero LT 25) then begin

		print," ...averaging..."
		mosaic_image = mosaic_image/(total_weight>1.e-37)

	  endif else begin

		wd = where( total_weight GT 0, npix )	;Locations of data
		print," ...averaging ", npix, " non-empty pixels..."
		mosaic_image(wd) = mosaic_image(wd) / total_weight(wd)
	   endelse

	if keyword_set( back ) then begin

		if N_elements( wd ) LE 0 then begin
			wd = where( total_weight GT 0, npix )
			print," ...averaged ", npix, " non-empty pixels..."
		   endif

		maxval = max( mosaic_image(wd), MIN=minval )
		print," MAX =",maxval,"    MIN =",minval
		range = maxval - minval
		zeroLevel = minval-.001*range
		print," setting background to:", zeroLevel
		mosaic_image(wz) = zeroLevel
	   endif

return, conv_vartype( mosaic_image, TYPE=type_out )
end