source: tags/release-0.9.2/Utils/wcsFunctions.cc @ 201

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

A swag of changes, centred around the addition of functionality to read
in an already saved FITS file containing the reconstructed array.

duchamp.hh -- const strings containing new keywords for saved FITS files.
Cubes/saveImage.cc -- improved methods for saving the FITS files, including

header keywords indicating origin.

Cubes/readRecon.cc -- new function to read in reconstructed array from FITS file.
Cubes/cubes.cc -- recon array now declared if recon file exists.
Cubes/cubes.hh -- prototype for readReconCube
Cubes/trimImage.cc -- allow for case of recon file existing.
Cubes/getImage.cc -- improved error/warning reporting to match new function.
param.hh -- new parameters to deal with reconFile
param.cc -- IO and default values for new parameters.
docs/Guide.tex -- Describing new option. Also added instal guide.
InputComplete? -- Added new parameters, and set values shown to be default values.
mainDuchamp.cc -- code related to new function.
ATrous/ReconSearch.cc -- will now not do reconstruction if recon array exists.

Fixed bug in output.

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