function ext_subimg3, index_in, data_in, index_out, $
	update_index=update_index, $
	ref_image=ref_image, $
	fri=fri0, nx=nx, ny=ny, center=center, $
	use_avg=use_avg, qprint=qprint, typout=typout, $
	margin=margin, trim=trim, $
	qdebug=qdebug
;
;+
;NAME:
;	ext_subimg3
;PURPOSE:
;	To extract a sub image out of an image.  The resolution of the image 
;	is not changed.
;METHOD:
;	There are two ways to extract.
;	     1. When asking for NxM, always return NxM no matter if the
;		input image does not fall within the coordinates requested
;		(or partially overlaps the coordinates requested)
;		This option is default.
;	     2. When asking for NxM, and only M1 lines of the input image
;		overlap the requested coordinates, then return an NxM1 sized
;		image.  This option is enabled by using /TRIM option
;SAMPLE CALLING SEQUENCE:
;	ext_subimg3, index, data, index2, data2, ref_image=pfi_index
;INPUT:
;	index_in - The index for the FFI image
;	data_in	- The data array of the FFI image
;OUTPUT:
;	index_out- The index for the output.  
;	data_out- The data array extracted from the FFI
;OPTIONAL KEYWORD INPUT:
;	fri	- The Full Resolution IDL coordinate of the lower left sub-pixel
;		  in the lower left of the observing region to extract. [x0,y0]
;		  NOTE: It is possible that the corner return will not be
;		  this exact value because HR and QR images exist only on
;		  every other (and every fourth) pixel.  The history record
;		  will reflect the actual portion extracted.
;	nx	- The number of full resolution pixels to extract in X
;	ny	- The number of full resolution pixels to extract in Y
;	center	- If set, then the "FRI" coordinates are for the center of the
;		  field of view to extract
;	ref_image- It is possible to pass a PFI index structure and have the
;		  corresponding portion extracted.
;	use_avg	- If set, and if the requested image falls outside of the
;		  input image, then set the portion outside to the average
;		  value of the input image (using the 50 lines closest to
;		  the output image).  Used for dark subtraction stuff.
;	margin	- Can be used in conjuntion with "ref_image".  A extra
;		  number of pixels will be extracted from around the image
;		  described by "ref_image".  The units are in the summed
;		  pixels.  It should be a scalar.
;	trim	- If set, then trim the image to only return the portion that
;		  overlaps the requested portion.
;RESTRICTIONS:
;	The routine assumes that it is working on raw FFI images, so that
;	they are complete in the E/W direction, but that the FR images
;	might be strips so they will not be complete in the N/S direction.
;HISTORY:
;	Written 6-Jun-93 by M.Morrison
;	20-Jun-93 (MDM) - Added more header information
;			- Inserted call to GT_CORNER
;	25-Aug-93 (MDM) - Whole series of major changes
;	26-Aug-93 (MDM) - More changes
;	 7-Oct-93 (MDM) - Patched bug causing array out of bounds
;	12-Oct-93 (MDM) - Patch for the case where trying to extract too far 
;			  to the left of the image
;	18-Oct-93 (MDM) - Added FRI=FIX(FRI) statement
;			- Corrected a roundoff problem (integer versus real)
;V1.02	24-Nov-93 (MDM) - Corrected an error.  It was returning a full 64x64
;			  image even when the overlap between the input image
;			  and the requested portion was something like 64x48
;			  lines.  The other 16 lines had zeros and was causing
;			  problems with SXT_PREP2.  Added /TRIM option.
;			- Also modified the SHAPE_SAV technique
;V1.03	29-Nov-93 (MDM) - Finished 24-Nov correction - it was causing a problem
;			  when the region requested was way off (x1>0, y1>0
;			  was coded)
;V1.04	16-Dec-93 (MDM) - Recent mods caused a certain condition not to be
;			  recognized.  Corrected it (added yy1 variable)
;-
;
progverno = 1.03*1000
qprint = keyword_set(qprint)
his_index, /enable
if (n_elements(margin) eq 0) then margin = 0
qtrim = keyword_set(trim)
;
;------------------------------------------------------------------------------
;		Definition of what to extract
;------------------------------------------------------------------------------
;
if (keyword_set(ref_image)) then begin
    fri = gt_corner(ref_image, /lower_left)	;fri(0) is now the lower left sub pixel of the summed pixels
    factor = 2.^gt_res(ref_image) / 2.^gt_res(index_in)		;pass in HR, reference FR means 1/2 as many pixels extracted
    if (n_elements(nx) eq 0) then nx = gt_shape(ref_image, /x)*factor + 2*margin
    if (n_elements(ny) eq 0) then ny = gt_shape(ref_image, /y)*factor + 2*margin
end
if (keyword_set(fri0)) then begin
    if (keyword_set(center)) then fri = fri0 - fix([nx,ny]/2) $
				else fri = fri0
end
fri = fix(fri)		;FRI should always be an integer value - MDM added 18-Nov-93
;
if (n_elements(fri) eq 0) then begin
    message, 'Corner must be specified either by FRI or REF_IMAGE keyword', /info
    return, 0
end
;
;
;------------------------------------------------------------------------------
;		Definition of input array
;------------------------------------------------------------------------------
;
siz = size(data_in)
typ_in = siz( siz(0)+1 )
nx_in = siz(1)
ny_in = siz(2)

;------------------------------------------------------------------------------
;		Extract
;------------------------------------------------------------------------------
;
res = fix(gt_res(index_in))	;made RES an integer - MDM 18-Nov-93
npix = 2^res	;1,2,4
xc = fri(0)	;x corner to be extracted (requested corner of the output array)
yc = fri(1)	;y corner to be extracted
;
ix0 = gt_corner(index_in, /x, /lower)		; corner address of input image
iy0 = gt_corner(index_in, /y, /lower)
;;nx_in = gt_shape(index_in, /x)		; MDM removed 24-Nov-93 - Use the data array size instead - avoid PFI/OR confusion
;;ny_in = gt_shape(index_in, /y)
ox0 = 0		;insertion coordinates into the output array
oy0 = 0
;
x0 = fix(xc-ix0)/npix - margin		;pixel location in the input array (in summed pixels) which should be extracted
y0 = fix(yc-iy0)/npix - margin
;
if (x0 lt 0) then begin		;MDM added 12-Oct-93
   ox0 = fix(abs(x0)+0.5)		;for the case where trying to extract too far to the left of the image
   x0 = 0
end
if ((y0 lt 0) and qtrim) then begin	;MDM added 24-Nov-93
    oy0 = fix(abs(y0))
    y0 = 0
end
x1 = (x0 + nx-1 - ox0)<(nx_in-1)>0
yy1 = (y0 + ny-1 - oy0)			;MDM added 16-Dec-93
y1 = yy1<(ny_in-1)>0
;
if (keyword_set(qdebug)) then begin
    print, 'Input image LL corner: ', ix0, iy0
    print, 'Requested LL corner:   ', xc, yc
    print, 'x0,y0,x1,y1:           ', x0,y0,x1,y1
    print, 'ox0,oy0:               ', ox0,oy0
end
;
;------------------------------------------------------------------------------
;		Definition of output array
;------------------------------------------------------------------------------
;
if (n_elements(typout) eq 0) then typout = typ_in
if (not qtrim) then data_out = make_array(nx, ny, type=typout)
qfloat = typout ge 4
;
;------------------------------------------------------------------------------
;		Check for output image falling outside of input image
;------------------------------------------------------------------------------

if ((y0 lt 0) or (yy1 ge ny_in) and (not qtrim)) then begin
    y00 = 0
    if (y0 lt 0) then begin
	y00 = abs(y0)  < (ny)
	if (keyword_set(use_avg)) then begin
	    temp = data_in(*, 0:49)
	    signal = total(temp)/n_elements(temp)
	    if (not qfloat) then signal = fix(signal + 0.5)
	    data_out(*,0:y00-1) = signal
	    if (qprint) then message, 'Setting lines 0 to ' + strtrim(y00-1,2)+ ' equal to '+ strtrim(signal,2)+ ' DN', /info
	end
    end

    if ((y0 lt ny_in) and (yy1 ge 0)) then data_out(ox0,y00) = data_in(x0:x1, y0>0:yy1<(ny_in-1))

    if (yy1 ge ny_in) then begin
	n_over = yy1-(ny_in-1) + ((y0-ny_in)>0)
	y11 = (ny - n_over)>0		;>0 added 12-May-93
	if (keyword_set(use_avg)) then begin
	    temp = data_in(*, ny_in-50:ny_in-1)
	    signal = total(temp)/n_elements(temp)
	    if (not qfloat) then signal = fix(signal + 0.5)
	    data_out(*,y11:ny-1) = signal
	    if (qprint) then message, 'Setting lines ' + strtrim(y11,2) + ' to ' + $
								strtrim(ny-1,2)+ ' equal to '+ strtrim(signal,2)+ ' DN', /info
	end
    end
end else begin
    if (qtrim) then data_out = data_in(x0:x1, y0:y1) $
		else data_out(ox0,oy0) = data_in(x0:x1, y0:y1)
end

index_out = index_in
corner = gt_corner(index_in)	; corner address of input image - center of summed pixel
if (qtrim) then corner = corner + [x0,y0]*npix	$	; shift corner by x0,y0 which is the location where the extraction occurred
	else corner = corner + [x0-ox0,y0-oy0]*npix
;;shape = [nx, ny]
shape = [x1-x0+1, y1-y0+1]	;MDM modified 24-Nov-93
his_index, index_out, 0, 'corner_ext', corner
his_index, index_out, 0, 'corner_sav', corner
index_out.sxt.shape_sav = shape
;
if (keyword_set(qdebug)) then print, 'History Corner: ', corner
if (keyword_set(update_index)) then index_in = index_out
if (keyword_set(qstop)) then stop
;
return, data_out
end