pro ssw_reg, iwcs, idata, owcs, odata, wcs_ref=wcs_ref, $
  scale_ref=scale_ref, roll_ref=roll_ref, crpix1_ref, crpix2_ref, $
  x0_ref=x0_ref, y0_ref=y0_ref, $
  interp=interp, cubic=cubic, x_off=x_off, y_off=y_off, $
  normalize=normalize, $
  no_remap=no_remap, $
  _extra=_extra, qstop=qstop

;+
; NAME:
;   ssw_reg
; PURPOSE:
;   Register any image (rotate, translate, scale) to arbitrary reference values of
;   scale, roll, and image center x and y, given a structure with this minimal tag set:
;     [naxis1, naxis2, crpix1, crpix2, cdelt1, cdelt2, crota2]
; CATEGORY:
;   Image registration
; CALL:
;
; OUTPUTS:
;   owcs - The updated index structure of the input image
;   odata - Registered output image.
; KEYWORDS:
; TODO:
; HISTORY:
;   2012-11-15 - GLS - Fixed errors in definition of CRPIX1,2 in OINDEX.
;-

; Define prognam, progver variables
prognam = 'SSW_REG.PRO'

;progver = 'V1.0' ; 2012-10-17 (GLS)
;progver = 'V1.1' ; 2012-11-15 (GLS)
progver = 'V1.2' ; 2013-05-28 (GLS)
progver = 'V1.3' ; 2013-06-27 (GLS)

; Define defaults for ROT interpolation and fill parameters:
if ( (not exist(interp)) and (not exist(cubic)) ) then cubic = -0.5
if not exist(missing) then missing = min(idata)

if n_params() eq 0 then begin
  print, ' Program version number = ' + strtrim(progver,2)
  return
endif

if not exist(roll_ref) then roll_ref = wcs_ref.crota2
if not exist(scale_ref) then scale_ref = wcs_ref.cdelt1

img1_cen_pix = wcs_pix2pix(wcs_ref, iwcs)
if exist(crpix1_ref) then x0 = crpix1_ref-1.0 else x0 = img1_cen_pix[0]
if exist(crpix2_ref) then y0 = crpix2_ref-1.0 else y0 = img1_cen_pix[1]

; Define deltas for rotation, translation, and scaling by subtracting the reference
; parameters from the corresponding image parameters:
delta_roll = roll_ref - iwcs.crota2
scale_fac = iwcs.cdelt1/scale_ref

if not keyword_set(no_remap) then begin

; Call ROT.PRO to re-map image:
;  odata = rot(idata, delta_roll, scale_fac, crpix1_ref-1.0, crpix2_ref-1.0, $
;              interp=interp, cubic=cubic, pivot=pivot, missing=missing)
   odata = rot(idata, delta_roll, scale_fac, x0, y0, $
               interp=interp, cubic=cubic, pivot=pivot, missing=missing)

endif else begin

; Allow 'no_remap' for testing purposes:
   odata = idata

endelse

; Set coordinate tags equal to reference values:
owcs = iwcs
owcs.crota2 = roll_ref
owcs.cdelt1 = scale_ref
owcs.cdelt2 = scale_ref
; Update CRPIX1,2, including offsets for differences in sizes of
; input and reference image arrays:
owcs.crpix1 = (iwcs.naxis1/2.-wcs_ref.naxis1/2.) + wcs_ref.crpix1
owcs.crpix2 = (iwcs.naxis2/2.-wcs_ref.naxis2/2.) + wcs_ref.crpix2

if keyword_set(normalize) and tag_exist(iwcs, 'exptime') then begin
  odata = odata/iwcs.exptime
  owcs.exptime = 1.0
endif

if tag_exist(owcs, 'lvl_num') then owcs.lvl_num = 1.5

; Update history record:

buff_history = ''

interpolation_algorithm = 'ROT called with nearest neighbor sampling'
if exist(interp) then $
  if interp ne 0 then $
    interpolation_algorithm = 'ROT called with linear interpolation'
if exist(cubic) then $
  if cubic ne 0 then $
    interpolation_algorithm = 'ROT called with cubic interpolation: cubic = ' + strtrim(cubic,2)

buff_history = [buff_history, interpolation_algorithm]

if ~tag_exist(wcs_ref, 'pnt_ref') then begin
  if tag_exist(wcs_ref,'fsn') then $
    buff_history = [buff_history, 'Image registered to SDO image with FSN = ' + $
                                   strtrim(wcs_ref.fsn,2)]
  if tag_exist(wcs_ref,'t_obs') then $
    buff_history = [buff_history, 'Image registered to SDO image with T_OBS = ' + $
                                   strtrim(wcs_ref.t_obs,2)]
endif else begin
  buff_history = [buff_history, 'Image registered to sun center, roll 0, and ' + $
                                'platescale = 0.6 arcsec per pixel.']
endelse

update_history, owcs, version=progver, caller=prognam
update_history, owcs, buff_history

if keyword_set(qstop) then stop,' Stopping on request.'

end