1
2
3
4
5
6 package org.astrogrid.datacenter.fits;
7
8
9 /*** Calculates World COordinates from the given FITS header.
10 *
11 * NB THIS IS BLATENTLY STOLEN FROM Pierre Grosbol TO TEST WITH, THIS NEEDS
12 * TO BE SORTED BEFORE RELEASE!
13 *
14 * @version $Revision: 1.1 $ $Date: 2004/09/28 15:02:13 $
15 * @author P.Grosbol, DMD/ESO, <pgrosbol@eso.org>
16 */
17 public class FitsWCS {
18
19 public final static int LIN = 0;
20 public final static int TAN = 1;
21 public final static int ARC = 2;
22
23 private final static int MPS = 20;
24
25
26 protected int nax = 0;
27 protected int[] cproj;
28 protected double[] crpix;
29 protected double[] crval;
30 protected double[] cdelt;
31 protected double[] crota;
32 protected String[] ctype;
33 protected double[][] cdMatrix;
34 protected double[][] pcMatrix;
35 protected boolean hasPcMatrix = false;
36 protected boolean hasCdMatrix = false;
37
38 protected double[] amdx = new double[MPS];
39 protected double[] amdy = new double[MPS];
40
41 /*** Default constructor for FitsWCS class.
42 *
43 public FitsWCS() {
44 }
45
46 /** Constructor for FitsData class given a FITS header with
47 * associated data unit as a file.
48 *
49 * @param header FitsHeader object with the image header
50 */
51 public FitsWCS(FitsHeader header) {
52 this(header, ' ');
53 }
54
55 /*** Constructor for FitsData class given a FITS header with
56 * associated data unit as a file.
57 *
58 * @param header FitsHeader object with the image header
59 * @param ver version of WCS i.e. ' ' or 'A'..'Z'
60 */
61 public FitsWCS(FitsHeader header, char ver) {
62 setHeader(header, ver);
63 }
64
65
66 /*** Define FITS header for FitsWCS object.
67 *
68 * @param header FitsHeader object with the image header
69 * @param ver version of WCS i.e. ' ' or 'A'..'Z'
70 */
71 private void setHeader(FitsHeader header, char ver) {
72
73 FitsKeyword kw = header.get("NAXIS");
74 nax = (kw == null) ? 0 : kw.toInt();
75 init(nax);
76
77 String sver = String.valueOf(ver).toUpperCase().trim();
78 for (int j=1; j<=nax; j++) {
79 kw = header.get("CRPIX"+j+sver);
80 crpix[j-1] = (kw == null) ? 0.0 : kw.toReal();
81
82 kw = header.get("CRVAL"+j+sver);
83 crval[j-1] = (kw == null) ? 0.0 : kw.toReal();
84
85 kw = header.get("CDELT"+j+sver);
86 cdelt[j-1] = (kw == null) ? 1.0 : kw.toReal();
87 cdMatrix[j-1][j-1] = cdelt[j-1];
88
89 kw = header.get("CTYPE"+j+sver);
90 ctype[j-1] = (kw == null) ? " " : kw.getValue();
91 if (7<ctype[j-1].length()) {
92 String wctype = ctype[j-1].substring(5, 7);
93 if (wctype.equals("TAN")) {
94 cproj[j-1] = TAN;
95 } else if (wctype.equals("ARC")) {
96 cproj[j-1] = ARC;
97 } else {
98 cproj[j-1] = LIN;
99 }
100 } else {
101 cproj[j-1] = LIN;
102 }
103
104 for (int i=1; i<=nax; i++) {
105 kw = header.get("CD"+j+"_"+j+sver);
106 if (kw != null) {
107 cdMatrix[j-1][i-1] = kw.toReal();
108 hasCdMatrix = true;
109 }
110 }
111
112 for (int i=1; i<=nax; i++) {
113 kw = header.get("PC"+j+"_"+i+sver);
114 if (kw != null) {
115 pcMatrix[j-1][i-1] = kw.toReal();
116 hasPcMatrix = true;
117 }
118 }
119 }
120
121 for (int j=1; j<MPS; j++) {
122 kw = header.get("AMDX"+j);
123 amdx[j-1] = (kw == null) ? 0.0 : kw.toReal();
124 kw = header.get("AMDY"+j);
125 amdy[j-1] = (kw == null) ? 0.0 : kw.toReal();
126 }
127 }
128
129 /*** Initiate internal WCS data structures.
130 *
131 * @param nax No. of axies of data array
132 */
133 private void init(int nax) {
134 cproj = new int[nax];
135 crpix = new double[nax];
136 crval = new double[nax];
137 cdelt = new double[nax];
138 ctype = new String[nax];
139 cdMatrix = new double[nax][nax];
140 pcMatrix = new double[nax][nax];
141 ctype = new String[nax];
142
143 for (int n=0; n<nax; n++) {
144 cproj[n] = LIN;
145 crpix[n] = 0.0;
146 crval[n] = 0.0;
147 cdelt[n] = 1.0;
148 ctype[n] = " ";
149 for (int i=0; i<nax; i++) {
150 cdMatrix[n][i] = (i==n) ? 1.0 : 0.0;
151 pcMatrix[n][i] = (i==n) ? 1.0 : 0.0;
152 }
153 }
154 }
155
156 /*** Compute World Coordinates from pixel coordinates.
157 *
158 * @param pix Array with pixel coordinates
159 */
160 public double[] toWCS(double[] p) {
161 double[] w;
162 double[] q;
163
164 for (int j=0; j<nax; j++) p[j] -= crpix[j];
165 if (hasPcMatrix) {
166 q = matrixMult(pcMatrix, p);
167 for (int j=0; j<nax; j++) q[j] *= cdelt[j];
168 } else if (hasCdMatrix) {
169 q = matrixMult(cdMatrix, p);
170 } else {
171 q = p;
172 for (int j=0; j<nax; j++) q[j] *= cdelt[j];
173 }
174
175 switch (cproj[0]) {
176 case TAN:
177 case LIN:
178 for (int j=0; j<nax; j++) q[j] += crval[j];
179 }
180
181 return q;
182 }
183 /*** Compute pixel coordinates from a set of World Coordinates.
184 *
185 * @param wcs Array with World Coordinates
186 */
187 public double[] toPixel(double[] wcs) {
188 double[] pix = new double[nax];
189 double[] wc = new double[nax];
190
191 for (int n=0; n<nax; n++) {
192 pix[n] = (wcs[n]-crval[n])/cdelt[n] + crpix[n];
193 }
194 return pix;
195 }
196
197 private double[] matrixMult(double[][] mtx, double[] vec) {
198 int nv = vec.length;
199 double[] res = new double[nv];
200
201 for (int j=0; j<nv; j++) {
202 res[j] = 0;
203 for (int i=0; i<nv; i++) {
204 res[j] += mtx[j][i] * vec[i];
205 }
206 }
207 return res;
208 }
209 }
210
211
212
213
214