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