001 /* @(#)FitsWCS.java $Revision: 1.3 $ $Date: 2003/04/08 14:54:31 $
002 *
003 * Copyright (C) 2003 European Southern Observatory
004 * License: GNU General Public License version 2 or later
005 */
006 package org.eso.fits;
007
008 import java.io.*;
009 import java.text.*;
010
011 /** FitsData class represents a FITS data unit
012 *
013 * @version $Revision: 1.3 $ $Date: 2003/04/08 14:54:31 $
014 * @author P.Grosbol, DMD/ESO, <pgrosbol@eso.org>
015 */
016 public class FitsWCS {
017
018 public final static int LIN = 0;
019 public final static int TAN = 1;
020 public final static int ARC = 2;
021
022 private final static int MPS = 20;
023
024 protected int type;
025 protected int nax = 0;
026 protected int[] cproj;
027 protected double[] crpix;
028 protected double[] crval;
029 protected double[] cdelt;
030 protected double[] crota;
031 protected String[] ctype;
032 protected double[][] cdMatrix;
033 protected double[][] pcMatrix;
034 protected boolean hasPcMatrix = false;
035 protected boolean hasCdMatrix = false;
036
037 protected double[] amdx = new double[MPS];
038 protected double[] amdy = new double[MPS];
039
040 /** Default constructor for FitsWCS class.
041 */
042 public FitsWCS() {
043 type = Fits.UNKNOWN;
044 }
045
046 /** Constructor for FitsData class given a FITS header with
047 * associated data unit as a file.
048 *
049 * @param header FitsHeader object with the image header
050 */
051 public FitsWCS(FitsHeader header) {
052 this();
053 setHeader(header, ' ');
054 }
055
056 /** Constructor for FitsData class given a FITS header with
057 * associated data unit as a file.
058 *
059 * @param header FitsHeader object with the image header
060 * @param ver version of WCS i.e. ' ' or 'A'..'Z'
061 */
062 public FitsWCS(FitsHeader header, char ver) {
063 this();
064 setHeader(header, ver);
065 }
066
067 /** Constructor for FitsData class given a No. of axies of data array.
068 *
069 * @param nax No. of axies of data array
070 */
071 public FitsWCS(int nax) {
072 this();
073 init(nax);
074 }
075
076 /** Define FITS header for FitsWCS object.
077 *
078 * @param header FitsHeader object with the image header
079 * @param ver version of WCS i.e. ' ' or 'A'..'Z'
080 */
081 public void setHeader(FitsHeader header, char ver) {
082 type = header.getType();
083
084 FitsKeyword kw = header.getKeyword("NAXIS");
085 nax = (kw == null) ? 0 : kw.getInt();
086 init(nax);
087
088 String sver = String.valueOf(ver).toUpperCase().trim();
089 for (int j=1; j<=nax; j++) {
090 kw = header.getKeyword("CRPIX"+j+sver);
091 crpix[j-1] = (kw == null) ? 0.0 : kw.getReal();
092
093 kw = header.getKeyword("CRVAL"+j+sver);
094 crval[j-1] = (kw == null) ? 0.0 : kw.getReal();
095
096 kw = header.getKeyword("CDELT"+j+sver);
097 cdelt[j-1] = (kw == null) ? 1.0 : kw.getReal();
098 cdMatrix[j-1][j-1] = cdelt[j-1];
099
100 kw = header.getKeyword("CTYPE"+j+sver);
101 ctype[j-1] = (kw == null) ? " " : kw.getString();
102 if (7<ctype[j-1].length()) {
103 String wctype = ctype[j-1].substring(5, 7);
104 if (wctype.equals("TAN")) {
105 cproj[j-1] = TAN;
106 } else if (wctype.equals("ARC")) {
107 cproj[j-1] = ARC;
108 } else {
109 cproj[j-1] = LIN;
110 }
111 } else {
112 cproj[j-1] = LIN;
113 }
114
115 for (int i=1; i<=nax; i++) {
116 kw = header.getKeyword("CD"+j+"_"+j+sver);
117 if (kw != null) {
118 cdMatrix[j-1][i-1] = kw.getReal();
119 hasCdMatrix = true;
120 }
121 }
122
123 for (int i=1; i<=nax; i++) {
124 kw = header.getKeyword("PC"+j+"_"+i+sver);
125 if (kw != null) {
126 pcMatrix[j-1][i-1] = kw.getReal();
127 hasPcMatrix = true;
128 }
129 }
130 }
131
132 for (int j=1; j<MPS; j++) {
133 kw = header.getKeyword("AMDX"+j);
134 amdx[j-1] = (kw == null) ? 0.0 : kw.getReal();
135 kw = header.getKeyword("AMDY"+j);
136 amdy[j-1] = (kw == null) ? 0.0 : kw.getReal();
137 }
138 }
139
140 /** Initiate internal WCS data structures.
141 *
142 * @param nax No. of axies of data array
143 */
144 private void init(int nax) {
145 cproj = new int[nax];
146 crpix = new double[nax];
147 crval = new double[nax];
148 cdelt = new double[nax];
149 ctype = new String[nax];
150 cdMatrix = new double[nax][nax];
151 pcMatrix = new double[nax][nax];
152 ctype = new String[nax];
153
154 for (int n=0; n<nax; n++) {
155 cproj[n] = LIN;
156 crpix[n] = 0.0;
157 crval[n] = 0.0;
158 cdelt[n] = 1.0;
159 ctype[n] = " ";
160 for (int i=0; i<nax; i++) {
161 cdMatrix[n][i] = (i==n) ? 1.0 : 0.0;
162 pcMatrix[n][i] = (i==n) ? 1.0 : 0.0;
163 }
164 }
165 }
166
167 /** Compute World Coordinates from pixel coordinates.
168 *
169 * @param p Array with pixel coordinates
170 */
171 public double[] toWCS(double[] p) {
172 double[] w;
173 double[] q;
174
175 for (int j=0; j<nax; j++) p[j] -= crpix[j];
176 if (hasPcMatrix) {
177 q = matrixMult(pcMatrix, p);
178 for (int j=0; j<nax; j++) q[j] *= cdelt[j];
179 } else if (hasCdMatrix) {
180 q = matrixMult(cdMatrix, p);
181 } else {
182 q = p;
183 for (int j=0; j<nax; j++) q[j] *= cdelt[j];
184 }
185
186 switch (cproj[0]) {
187 case TAN:
188 case LIN:
189 for (int j=0; j<nax; j++) q[j] += crval[j];
190 }
191
192 return q;
193 }
194 /** Compute pixel coordinates from a set of World Coordinates.
195 *
196 * @param wcs Array with World Coordinates
197 */
198 public double[] toPixel(double[] wcs) {
199 double[] pix = new double[nax];
200 double[] wc = new double[nax];
201
202 for (int n=0; n<nax; n++) {
203 pix[n] = (wcs[n]-crval[n])/cdelt[n] + crpix[n];
204 }
205 return pix;
206 }
207
208 private double[] matrixMult(double[][] mtx, double[] vec) {
209 int nv = vec.length;
210 double[] res = new double[nv];
211
212 for (int j=0; j<nv; j++) {
213 res[j] = 0;
214 for (int i=0; i<nv; i++) {
215 res[j] += mtx[j][i] * vec[i];
216 }
217 }
218 return res;
219 }
220 }
221
222
223
224
225