source: trunk/Utils/wcsFunctions.cc @ 22

Last change on this file since 22 was 22, checked in by Matthew Whiting, 18 years ago

Corrected the pix--wcs transformations to account for an offset of 1 pixel,
due to my arrays being indexed to 0, and FITS arrays indexed to 1.
New wcsToPix/pixToWCS functions written to handle multiple pixels, and
code added to make use of them in plotting.cc and detection.cc.
Code added to make use of setVel_kms in same files.
File size: 6.3 KB
Line 
1#include <iostream>
2#include <wcs.h>
3#include <Utils/utils.hh>
4
5double pixelToVelocity(wcsprm *wcs, double &x, double &y, double &z)
6{
7  /**
8   *  pixelToVelocity(wcsprm *wcs, double &x, double &y, double &z)
9   *   Uses wcs to convert the three-dimensional pixel position (x,y,z)
10   *   to world coordinates.
11   *   Returns the velocity.
12   */
13
14  int    *stat   = new int[1];
15  double *pixcrd = new double[3];
16  double *imgcrd = new double[3];
17  double *world  = new double[3];
18  double *phi    = new double[1];
19  double *theta  = new double[1];
20  pixcrd[0] = x;
21  pixcrd[1] = y;
22  pixcrd[2] = z;
23  if(int flag=wcsp2s(wcs, 1, 3, pixcrd, imgcrd, phi, theta, world, stat)>0){
24    std::cerr<<"WCS Error in pixelToVelocity(): Code = "<<flag<<": "<<wcs_errmsg[flag]<<std::endl;
25    std::cerr<<"  stat value is "<<stat[0]<<std::endl;
26  }
27
28  double vel = setVel_kms(wcs, world[2]);
29
30  delete [] stat;
31  delete [] pixcrd;
32  delete [] imgcrd;
33  delete [] world;
34  delete [] phi;
35  delete [] theta;
36 
37  return vel;
38}
39 
40
41int pixToWCSSingle(wcsprm *wcs, const double *pix, double *world)
42{
43  /**
44   *  pixToWCSSingle(wcsprm *wcs, const double *pix, double *world)
45   *   Uses wcs to convert the three-dimensional pixel position referenced by pix
46   *   to world coordinates, which are placed in the array world[].
47   *   Assumes that pix only has one point with an x,y,and z pixel positions.
48   *   Offsets these pixel values by 1 to account for the C arrays being indexed to 0.
49   */
50
51  int nelem=3,npts=1,flag;
52
53  double *newpix = new double[nelem*npts];
54  for(int i=0;i<nelem*npts;i++) newpix[i] = pix[i] + 1.; 
55  // correct from 0-indexed to 1-indexed pixel array
56
57  int    *stat   = new int[npts];
58  double *imgcrd = new double[nelem*npts];
59  double *phi    = new double[npts];
60  double *theta  = new double[npts];
61  if(flag=wcsp2s(wcs, npts, nelem, newpix, imgcrd, phi, theta, world, stat)>0){
62    std::cerr << "WCS Error in pixToWCSSingle(): Code = " <<flag<<": "<<wcs_errmsg[flag]<<std::endl;
63    std::cerr << "  stat value is " << stat[0] << std::endl;
64  }
65  delete [] stat;
66  delete [] imgcrd;
67  delete [] phi;
68  delete [] theta;
69  return flag;
70}
71
72int wcsToPixSingle(wcsprm *wcs, const double *world, double *pix)
73{
74  /**
75   *  wcsToPixSingle(wcsprm *wcs, const double *world, double *pix)
76   *   Uses wcs to convert the three-dimensional world coordinate position referenced by world
77   *   to pixel coordinates, which are placed in the array pix[].
78   *   Assumes that world only has one point (three values eg. RA, Dec, Velocity)
79   *   Offsets the pixel values by 1 to account for the C arrays being indexed to 0.
80   */
81
82  int nelem=3,npts=1,flag;
83  int    *stat   = new int[npts];
84  double *imgcrd = new double[nelem*npts];
85  double *phi    = new double[npts];
86  double *theta  = new double[npts];
87  if(flag=wcss2p(wcs, npts, nelem, world, phi, theta, imgcrd, pix, stat)>0){
88    std::cerr << "WCS Error in wcsToPixSingle(): Code = " <<flag<<": "<<wcs_errmsg[flag]<<std::endl;
89    std::cerr << "  stat value is " << stat[0] << std::endl;
90  }
91
92  for(int i=0;i<nelem;i++) pix[i] -= 1.;  // correct from 1-indexed to 0-indexed pixel array
93
94  delete [] stat;
95  delete [] imgcrd;
96  delete [] phi;
97  delete [] theta;
98  return flag;
99}
100
101int pixToWCSMulti(wcsprm *wcs, const double *pix, double *world, const int npts)
102{
103  /**
104   *  pixToWCSSingle(wcsprm *wcs, const double *pix, double *world)
105   *   Uses wcs to convert the three-dimensional pixel positions referenced by pix
106   *   to world coordinates, which are placed in the array world[].
107   *   pix is assumed to hold the positions of npts points.
108   *   Offsets these pixel values by 1 to account for the C arrays being indexed to 0.
109   */
110
111  int nelem=3,flag;
112
113  double *newpix = new double[nelem*npts];
114  for(int i=0;i<nelem*npts;i++) newpix[i] = pix[i] + 1.; 
115  // correct from 0-indexed to 1-indexed pixel array
116
117  int    *stat   = new int[npts];
118  double *imgcrd = new double[nelem*npts];
119  double *phi    = new double[npts];
120  double *theta  = new double[npts];
121  if(flag=wcsp2s(wcs, npts, nelem, newpix, imgcrd, phi, theta, world, stat)>0){
122    std::cerr << "WCS Error in pixToWCSSingle(): Code = " <<flag<<": "<<wcs_errmsg[flag]<<std::endl;
123    std::cerr << "  stat value is " << stat[0] << std::endl;
124  }
125  delete [] stat;
126  delete [] imgcrd;
127  delete [] phi;
128  delete [] theta;
129  return flag;
130}
131
132int wcsToPixMulti(wcsprm *wcs, const double *world, double *pix, const int npts)
133{
134  /**
135   *  wcsToPixSingle(wcsprm *wcs, const double *world, double *pix)
136   *   Uses wcs to convert the three-dimensional world coordinate position referenced by world
137   *   to pixel coordinates, which are placed in the array pix[].
138   *   world is assumed to hold the positions of npts points.
139   *   Offsets the pixel values by 1 to account for the C arrays being indexed to 0.
140   */
141  int nelem=3,flag;
142  int    *stat   = new int[npts];
143  double *imgcrd = new double[nelem*npts];
144  double *phi    = new double[npts];
145  double *theta  = new double[npts];
146  if(flag=wcss2p(wcs, npts, nelem, world, phi, theta, imgcrd, pix, stat)>0){
147    std::cerr << "WCS Error in wcsToPixSingle(): Code = " <<flag<<": "<<wcs_errmsg[flag]<<std::endl;
148    std::cerr << "  stat value is " << stat[0] << std::endl;
149  }
150
151  for(int i=0;i<nelem;i++) pix[i] -= 1.;  // correct from 1-indexed to 0-indexed pixel array
152
153  delete [] stat;
154  delete [] imgcrd;
155  delete [] phi;
156  delete [] theta;
157  return flag;
158}
159
160float setVel_kms(wcsprm *wcs, const double coord)
161{
162  /**
163   *  setVel_kms(wcsprm *wcs, const double coord)
164   *   Convert the wcs coordinate given by coord to a velocity in km/s.
165   *   Does this by checking the ztype parameter in wcs to see if it is FREQ or otherwise,
166   *   and converts accordingly.
167   */
168
169  string ztype = wcs->ctype[2];
170  if(ztype!="FREQ") {
171    return coord/1000.;
172  }
173  else{
174    return C_kms * (wcs->restfrq*wcs->restfrq - coord*coord) / (wcs->restfrq*wcs->restfrq + coord*coord);
175  }
176}
177
178double velToCoord(wcsprm *wcs, const float velocity)
179{
180  /**
181   *  velToCoord(wcsprm *wcs, const double coord)
182   *   Convert the velocity given to the appropriate world coordinate for wcs.
183   *   Does this by checking the ztype parameter in wcs to see if it is FREQ or otherwise,
184   *   and converts accordingly.
185   */
186
187  string ztype = wcs->ctype[2];
188  if(ztype!="FREQ")
189    return velocity * 1000.;
190  else
191    return wcs->restfrq * (1. - velocity/C_kms) / (1. + velocity/C_kms);
192 
193}
Note: See TracBrowser for help on using the repository browser.