function align_ar,datain,index,ref_image=ref_image,korrel=korrel_flag, $
ref_sunc=refsc, ref_xy=refxy,ref_ll=refll,ref_day=ref_day, $
missing=missing, ar_size=ar_size, center_lat=clat, $
center_lon=clon, center_date=cdate
ARC_SEC_PER_PIX = gt_pix_size() ; Used to convert to/from arseconds
nimage = n_elements(datain(0,0,*))
nx = n_elements(datain(*,0,0))
ny = n_elements(datain(0,*,0))
offset = fltarr(2,nimage)
if n_elements(ref_image) NE 1 then ref_image=nimage/2
ref_image = fix(ref_image)
if ref_image GT nimage-1 OR ref_image LT 0 then begin
message,'align_ar: ref_image is out of range'
endif
rd_pnt,index(0),index(0),pdata,flag=flag
if flag then message,'Could not find PNT file'
rd_fem,index(0),index(0),femdata,/near
for image=0,nimage-1 do begin
if int2secarr(index(image),pdata(0)) LT 0 OR $
int2secarr(index(image),pdata(n_elements(pdata)-1)) GT 0 $
then begin
; Read pointing information
rd_pnt,index(image),index(image),pdata,flag=flag
if flag then message,'Could not read PNT file'
endif
if int2secarr(index(image),femdata(0)) LT 0 OR $
int2secarr(index(image),femdata(n_elements(femdata)-1)) GT 0 $
then begin
; Read fem information
rd_fem,index(image),index(image),femdata,/near
endif
ii=tim2dset(pdata,index(image)) ; find the closest time to required time
; sc = gt_hxa(pdata(ii),/sxtpix,fm=fm) ; Get the sun center position
; if fm(0) then message,/info,'WARNING: HXA fiducial mark'
; sc = hxa_suncenter(pdata(ii),code=hxacode,fem=femdata) ; Get the sun center position
sc = get_suncenter(pdata,index=index(image),code=hxacode,fem=femdata) ; Get the sun center position
sc=[sc(0,0),sc(1,0)]
if hxacode(0) LE 1 then message,'WARNING: HXA Data Problem '+string(hxacode(0)),/info
if image EQ 0 then begin
; Get the reference heliographic coordinates and time
if int2secarr(index(ref_image),pdata(0)) LT 0 OR $
int2secarr(index(ref_image),pdata(n_elements(pdata)-1)) GT 0 $
then begin
rd_pnt,index(ref_image),index(ref_image),pdata,flag=flag
if flag then message,'Could not read PNT file'
endif
if int2secarr(index(ref_image),femdata(0)) LT 0 OR $
int2secarr(index(ref_image),femdata(n_elements(femdata)-1)) GT 0 $
then begin
rd_fem,index(ref_image),index(ref_image),femdata,/near
endif
ii=tim2dset(pdata,index(ref_image))
; refsc = gt_hxa(pdata(ii),/sxtpix) ; Get the sun center position
; refsc = hxa_suncenter(pdata(ii),code=hxacode,fem=femdata) ; Get the sun center position
refsc = get_suncenter(pdata,index=index(ref_image),code=hxacode,fem=femdata) ; Get the sun center position
refsc=[refsc(0,0),refsc(1,0)]
if hxacode(0) LE 1 then message,'WARNING: HXA Data Problem'+string(hxacode(0)),/info
sxt_fov_center, index(ref_image), refsc, xr,yr ; Rel. to SC
refxy = [xr,yr]
ll= arcmin2hel(xr(0)*ARC_SEC_PER_PIX/60.0,yr(0)*ARC_SEC_PER_PIX/60.0,$
date=fmt_tim(index(ref_image)))
latitude = ll(0) ; Set the initial lat/lon to center of first image
longitude = ll(1) ; Set the initial lat/lon to center of first image
ex2int,anytim2ex(index(ref_image)),msod,ds79
ref_day = ds79 + msod/8.64d7 ; AR reference time in decimal days
endif
; Use differential rotation to find the location of the active region
time = gt_time(index(image))
day = gt_day(index(image))
dday = (ref_day-(day+time/8.64d7)) ;difference in days from AR reference
sxt_fov_center, index(image), sc, xr, yr ; Rel. to SC
new_hel = arcmin2hel(xr(0)*ARC_SEC_PER_PIX/60.0, $
yr(0)*ARC_SEC_PER_PIX/60.0, $
date=fmt_tim(index(image)))
lat = new_hel(0)
lon = new_hel(1) + diff_rot(dday,lat)
; Where the center of the image *should* be (relative to sun center)
pix = hel2arcmin(lat,lon,date=fmt_tim(index(image)))
pix = pix*60.0/ARC_SEC_PER_PIX ; Full res pixels relative to sun center
; Compute the offset due to the solar rotation. This will be passed
; to mk_mosaic to shift all the images to the position (on the chip)
; of the reference image.
; In full resolution SXT pixels:
offset(0,image) = (pix(0)+refsc(0)) - (xr(0)+sc(0))
offset(1,image) = (pix(1)+refsc(1)) - (yr(0)+sc(1))
endfor
; Average value of the data array
if n_elements(missing) NE 1 then $
missing = total(datain)/n_elements(datain) $
else missing = missing(0)
; Make a full-res mosaic
data = mk_mosaic(datain,index,fix(offset),missing=missing(0), $
colstart=colstart,linstart=linstart)
nx = n_elements(data(*,0,0))
ny = n_elements(data(0,*,0))
; The reference coordinates passed out should refer to the center of the
; mosaic, not the center of the reference image
res = 2^gt_res(index(ref_image))
nxin = res*index(ref_image).sxt.shape_cmd(0)
nyin = res*index(ref_image).sxt.shape_cmd(1)
xmosaic_center = nx/2. - 0.5
ymosaic_center = ny/2. - 0.5
xin_center = nxin/2. - 0.5 + colstart(ref_image) + $
offset(0,ref_image) - fix(offset(0,ref_image))
yin_center = nyin/2. - 0.5 + linstart(ref_image) + $
offset(1,ref_image) - fix(offset(1,ref_image))
xcentoff = xmosaic_center - xin_center
ycentoff = ymosaic_center - yin_center
refxy(0) = refxy(0) + xcentoff
refxy(1) = refxy(1) + ycentoff
refll = arcmin2hel(refxy(0)*ARC_SEC_PER_PIX/60.0, $
refxy(1)*ARC_SEC_PER_PIX/60.0, $
date=fmt_tim(index(ref_image)))
d_xy = [0.,0.]
if n_elements(clat) GT 0 AND n_elements(clon) GT 0 then begin
; Adjust the PFI data cube so that the center is at the requested
; heliographic position
ex2int,anytim2ex(cdate),ms,ds
ddays = ref_day - (ds+ms/8.64d7)
latpfi = clat(0)
lonpfi = clon(0) + diff_rot(ddays,latpfi)
int2ex,long((ref_day-long(ref_day))*8.64d7),long(ref_day),pfidate
d_xy = -(hel2arcmin(latpfi,lonpfi,date=pfidate) $
- hel2arcmin(refll(0),refll(1),date=pfidate)) $
* 60.0/ARC_SEC_PER_PIX
if d_xy(0) NE 0.0 OR d_xy(1) NE 0.0 then refll=[latpfi,lonpfi]
endif
xoffset = 0.0
yoffset = 0.0
if n_elements(ar_size) EQ 2 then begin
; Extract a subarray to make the dimensions of data come out right
; Make sure that the center of the array remains the center of
; the array!
nxpfi = n_elements(data(*,0,0))
nypfi = n_elements(data(0,*,0))
if nxpfi NE ar_size(0) OR nypfi NE ar_size(1) then begin
data_size = size(data)
temp = make_array(ar_size(0),ar_size(1),n_elements(data(0,0,*)), $
type=data_size(n_elements(data_size)-2), $
value=missing)
xpfistart = (fix((ar_size(0)-nxpfi)/2.0)>0)
xpfiend = (fix((ar_size(0)+nxpfi)/2.0-1.0)<(ar_size(0)-1))
xpfirange = [(((nxpfi-ar_size(0))/2.0)>0), $
(((nxpfi+ar_size(0))/2.0-1.0)<(nxpfi-1))]
xoffset = xpfirange(0)-fix(xpfirange(0))
xpfirange = fix(xpfirange)
ypfistart = (fix((ar_size(1)-nypfi)/2.0)>0)
ypfiend = (fix((ar_size(1)+nypfi)/2.0-1.0)<(ar_size(1)-1))
ypfirange = [(((nypfi-ar_size(1))/2.0)>0), $
(((nypfi+ar_size(1))/2.0-1.0)<(nypfi-1))]
yoffset = ypfirange(0)-fix(ypfirange(0))
ypfirange = fix(ypfirange)
temp(xpfistart:xpfiend,ypfistart:ypfiend,*) = $
data(xpfirange(0):xpfirange(1), $
ypfirange(0):ypfirange(1),*)
data=temp & temp=0b
endif
endif
nx = n_elements(data(*,0,0))
ny = n_elements(data(0,*,0))
; The offset in the mosaic is only to the nearest pixel, so now shift
; the fractional part of the offset as well as any offsets from the
; lat/lon selection or the array size selection.
offset = offset-fix(offset)
for image=0,nimage-1 do begin
xshift = -offset(0,image)-xoffset-d_xy(0)
yshift = -offset(1,image)-yoffset-d_xy(1)
data(*,*,image) = poly_2d(data(*,*,image), $
[xshift,0,1,0],[yshift,1,0,0],1)
endfor
; Now the data is close, so do a correlation to get even closer.
; Each image is correlated with the previous image, which means a
; a cumulative shift from the first image must be computed later.
; The reason for doing it this way is that the individual shifts from
; the korrel program will be small and the korrel program can be run
; with a larger subimage and hence will run faster.
xshift = fltarr(nimage) ; Don't use the /nozero keyword here!!
yshift = fltarr(nimage)
if keyword_set(korrel_flag) then begin
for image=0,nimage-2 do begin
ii = (image+1+ref_image) MOD nimage ; First image
ii1 = (image+ref_image) MOD nimage ; Previous image
xrange = [nx/4,3*nx/4]
yrange = [ny/4,3*ny/4]
ref = reform(float(data(xrange(0):xrange(1),yrange(0):yrange(1),ii1)))
img = reform(float(data(xrange(0):xrange(1),yrange(0):yrange(1),ii)))
nxk = n_elements(ref(*,0)) ; This is the size of the mosaic'ed images
nyk = n_elements(ref(0,*))
MAX_SHIFT = 6 ; Should be an even positive integer
; Size of subimage in korrel:
sizes = [nxk-MAX_SHIFT,nyk-MAX_SHIFT]
; This centers the subimage:
offset = fix(0.5*([nxk,nyk]-sizes))
xy = korrel(img,ref,sizes,offset) ; Do the correlation
print,"korrel output: ",xy(0),xy(1)," image ",ii," of ",nimage
; Protect against a bad correlation
if max(abs(xy)) GT (MAX_SHIFT/2.0-1.0) then begin
xshift(ii) = 0.0
yshift(ii) = 0.0
endif $
else begin
xshift(ii) = xy(0,0) ; in full res pixels
yshift(ii) = xy(1,0) ; in full res pixels
endelse
endfor
; Now do the cumulative korrel shifts
xs = 0.0
ys = 0.0
for image=0,nimage-2 do begin
ii = (image+1+ref_image) MOD nimage
xs = xs + xshift(ii) ; Compute the cumulative shift
ys = ys + yshift(ii)
data(*,*,ii) = poly_2d(reform(data(*,*,ii)), $
[xs,0,1,0],[ys,1,0,0],1,nx,ny)
endfor
endif
return,data
end