-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathhrotate2.pro
212 lines (177 loc) · 7.31 KB
/
hrotate2.pro
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
pro hrotate2, oldim, oldhd, newim, newhd, direction,ERRMSG = errmsg
;+
; NAME:
; HROTATE
; PURPOSE:
; Apply the IDL ROTATE function and update astrometry in a FITS header
; EXPLANATION:
; Apply the intrinsic IDL ROTATE function to an image and update
; astrometry in the associated FITS header.
;
; CALLING SEQUENCE:
; HROTATE, oldim, oldhd, newim, newhd, direction
; or
; HROTATE, oldim, oldhd, direction
;
; INPUTS:
; OLDIM - the original image array
; OLDHD - the original FITS image header, string array
; DIRECTION - Scalar integer (0-7) specifying rotation direction,
; exactly as specified by the IDL ROTATE function.
;
; Direction Transpose? Rot. CCW X1 Y1
; ----------------------------------------
; 0 No None X0 Y0 (no change)
; 1 No 90 -Y0 X0
; 2 No 180 -X0 -Y0
; 3 No 270 Y0 -X0
; 4 Yes None Y0 X0
; 5 Yes 90 -X0 Y0
; 6 Yes 180 -Y0 -X0
; 7 Yes 270 X0 -Y0
;
; OPTIONAL OUTPUTS:
; NEWIM - the rotated image, with the same dimensions as Oldim
; NEWHD - header for newim containing updated astrometry info
; If output parameters are not supplied, the program
; will modify the input parameters OLDIM and OLDHD
; to contain the rotated image and updated header.
;
; OPTIONAL KEYWORD OUTPUT:
; ERRMSG - If this keyword is supplied, then any error mesasges will be
; returned to the user in this parameter rather than depending on
; on the MESSAGE routine in IDL. If no errors are encountered
; then a null string is returned.
; EXAMPLE:
; Rotate an image exactly 90 degrees counterclockwise and update the
; FITS image array and header.
;
; IDL> HROT, im, h, im_new, h_new, 1
;
; PROCEDURE:
; The image array is rotated using the ROTATE function.
; The CD (or CROTA) and CRPIX parameters, if present in the FITS header,
; are updated for the new rotation.
; History records are also added to the header
;
; RESTRICTIONS:
; Does not work Guide Star Survey (GSS) astrometry. Use GSSS_STDAST to
; first convert
; PROCEDURES USED:
; CHECK_FITS(), SXADDPAR, EXTAST
;
; MODIFICATION HISTORY:
; Written, Mar 1997 W. Landsman, Hughes STX
; Work for non-square images W. Landsman June 1998 Raytheon STX
; Converted to IDL V5.0 W. Landsman June 1998
; Fix for different plate scales, and CROTA2 defined, November 1998
; Added ERRMSG, Use double precision formatting, W. Landsman April 2000
; Consistent conversion between CROTA and CD matrix W. Landsman Oct 2000
; Correct update when CROTA keyword present W. Landsman June 2003
; Update CDELT for AIPS-style astrometry headers M. Perrin/WL Jul 2003
; Convert GSS astrometry to WCS W. Landsman November 2004
;-
On_error,2
npar = N_params()
if (npar NE 3) and (npar NE 5) then begin ;Check # of parameters
print,'Syntax - HROTATE, oldim, oldhd, newim, newhd, direction'
print,' or '
print,' HROTATE, oldim, oldhd, direction, {ERRMSG = ]'
return
endif
if npar EQ 3 then direction = newim
if N_elements(direction) NE 1 then message, $
'ERROR - Direction parameter must be an integer scalar (0-7)'
dirpar = direction mod 8
if dirpar LT 0 then dirpar = dirpar + 8
; Check that input header matches input image
save_err = arg_present(errmsg) ;Does user want error msgs returned?
; Check for valid 2-D image & header
; check_FITS, oldim, oldhd, dimen, /NOTYPE, ERRMSG = errmsg
; if errmsg NE '' then begin
; if not save_err then message,'ERROR - ' + errmsg,/CON
; return
; endif
dimen=[[(size(oldim))(1)],[(size(oldim))(2)]]
if N_elements(dimen) NE 2 then begin
errmsg = 'ERROR1 - Input image array must be 2-dimensional'
if not save_err then message,'ERROR - ' + errmsg,/CON
return
endif
if N_elements(dimen) NE 2 then message, $
'ERROR2 - Input image array must be 2-dimensional'
xsize = dimen[0] & ysize = dimen[1]
xc = (xsize-1)/2.
yc = (ysize-1)/2.
newhd = oldhd
if npar EQ 5 then newim = rotate(oldim, direction ) else $
oldim = rotate(oldim, direction )
if dirpar EQ 0 then return
label = 'HROTATE: ' + strmid(systime(),4,20)
sxaddhist, label + ' Image = ROTATE(Image,' + strtrim(direction,2) + ')',newhd
; Update astrometry info if it exists. If GSS astrometry is present, then
; convert it to standard WCS astrometry
extast, oldhd, astr, noparams
if strmid(astr.ctype[0],5,3) EQ 'GSS' then begin
gsss_stdast, newhd
extast, newhd, astr, noparams
endif
if noparams GE 0 then begin ;Astrometry parameters exist in header?
case dirpar of
0: rot_mat = [ [1,0], [0, 1] ]
1: rot_mat = [ [0, 1],[-1, 0] ]
2: rot_mat = [ [-1,0],[ 0,-1] ]
3: rot_mat = [ [0,-1], [1, 0] ]
4: rot_mat = [ [0, 1], [-1,0] ]
5: rot_mat = [ [-1,0], [0, -1] ]
6: rot_mat = [ [0,-1], [1, 0] ]
7: rot_mat = [ [1, 0], [0, 1] ]
else: message,$
'ERROR - Illegal value of direction parameter, must be between 0 and 7'
endcase
; For non-square images, check if X and Y axes have been flipped
crpix = astr.crpix
cd = astr.cd
cdelt = astr.cdelt
if cdelt[0] NE 1.0 then begin
cd[0,0] = cd[0,0]*cdelt[0] & cd[0,1] = cd[0,1]*cdelt[0]
cd[1,1] = cd[1,1]*cdelt[1] & cd[1,0] = cd[1,0]*cdelt[1]
endif
ncrpix = [xc,yc] + rot_mat#(crpix-1 -[xc,yc]) + 1
newcd = cd # transpose(rot_mat)
if (dirpar EQ 4) or (dirpar EQ 6) then begin
ncrpix[0] = xsize - ( ncrpix[0] - 1)
newcd[*,0] = -newcd[*,0]
endif
if (dirpar EQ 5) or (dirpar EQ 7) then begin
ncrpix[1] = ysize - (ncrpix[1] -1 )
newcd[*,1] = -newcd[*,1]
endif
if (xsize NE ysize) and (rot_mat[0,0] EQ 0) then begin
sxaddpar, newhd, 'NAXIS1', ysize
sxaddpar, newhd, 'NAXIS2', xsize
ncrpix[0] = ncrpix[0] - xc + yc
ncrpix[1] = ncrpix[1] - yc + xc
endif
sxaddpar, newhd, 'CRPIX1', ncrpix[0]
sxaddpar, newhd, 'CRPIX2', ncrpix[1]
if noparams EQ 3 then begin ;Transformation matrix format
sxaddpar, newhd, 'PC1_1', newcd[0,0]
sxaddpar, newhd, 'PC1_2', newcd[0,1]
sxaddpar, newhd, 'PC2_1', newcd[1,0]
sxaddpar, newhd, 'PC2_2', newcd[1,1]
endif else if noparams EQ 2 then begin
sxaddpar, newhd, 'CD1_1', newcd[0,0]
sxaddpar, newhd, 'CD1_2', newcd[0,1]
sxaddpar, newhd, 'CD2_1', newcd[1,0]
sxaddpar, newhd, 'CD2_2', newcd[1,1]
endif else begin ; noparams = 1. CROTA+CDELT type
crota = atan(-newcd[1,0], newcd[1,1] )*180.0/!DPI
if dirpar GE 4 then sxaddpar, newhd, 'CDELT1', -cdelt[0]
sxaddpar, newhd,'CROTA1', crota
sxaddpar, newhd,'CROTA2', crota
endelse
endif
if npar EQ 3 then oldhd = newhd ;update old image and header
return
end