source: tags/release-0.9/Utils/wcsFunctions.cc @ 813

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

Utils/wcsFunctions.cc Added a single-pixel-offset that had been missed in

rev.22.

Utils/position_related.cc Improved decToDMS to take into account different

axis types, and format the string accordingly.

Utils/utils.hh Corrected header call to match decToDMS change.
Detection/detection.cc Changed calls of decToDMS. Removed use of

Detection::lngRatio, as now unnecessary.

Detection/detection.hh Removed lngRatio and accessor functions from

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