;select subimage from an input image given xstart and ystart : pixels
;are symetrically removed from the input image according to these
;parameters

;astr is a structure containing the WCS parameters (from extast procedure)
pro select_rect_image, image, xstart, ystart, res, astr
nx_image = n_elements(image[*,0])
ny_image = n_elements(image[0,*])

;keep only pixels in center of image according to xystart
;size of output image
nx = nx_image -2*xstart
ny = ny_image -2*ystart

;shift to get to fist pixel
shift = ny_image*ystart+xstart
;1D array of subscripts given the size of the output image
use_x = shift+lindgen(nx)
;replicate it to get an array
use_xy = cmreplicate(use_x, [ny]) 

;create additionnal shifts (in x direction for each line)
shifts   = (2*xstart+nx)*lindgen(nx)
shifts_x = transpose(cmreplicate(shifts, [ny]))

use = shifts_x + use_xy
;print, shifts_x
res = dblarr(nx, ny)

res[*] = image[use[*]]

;update WCS parameters (pixels coordinates of reference point)
astr.crpix[0] = (nx + 1)/2.
astr.crpix[1] = (ny + 1)/2.

;stop
;return, res
end