[3] | 1 | #include <iostream> |
---|
[120] | 2 | #include <sstream> |
---|
[103] | 3 | #include <math.h> |
---|
[3] | 4 | #include <wcs.h> |
---|
[103] | 5 | #include <wcsunits.h> |
---|
[120] | 6 | #include <duchamp.hh> |
---|
[3] | 7 | #include <Utils/utils.hh> |
---|
| 8 | |
---|
| 9 | int pixToWCSSingle(wcsprm *wcs, const double *pix, double *world) |
---|
| 10 | { |
---|
[22] | 11 | /** |
---|
| 12 | * pixToWCSSingle(wcsprm *wcs, const double *pix, double *world) |
---|
[160] | 13 | * Uses wcs to convert the three-dimensional pixel position referenced |
---|
| 14 | * by pix to world coordinates, which are placed in the array world[]. |
---|
[22] | 15 | * Assumes that pix only has one point with an x,y,and z pixel positions. |
---|
[160] | 16 | * If there are any other axes present (eg. Stokes), these are set to the |
---|
| 17 | * first pixel in that axis. |
---|
| 18 | * Offsets these pixel values by 1 to account for the C arrays being |
---|
| 19 | * indexed to 0. |
---|
[22] | 20 | */ |
---|
| 21 | |
---|
[160] | 22 | int naxis=wcs->naxis,npts=1,flag; |
---|
[22] | 23 | |
---|
[160] | 24 | double *newpix = new double[naxis*npts]; |
---|
[22] | 25 | // correct from 0-indexed to 1-indexed pixel array |
---|
[160] | 26 | for(int i=0;i<naxis;i++) newpix[i] = 1.; |
---|
| 27 | newpix[wcs->lng] += pix[0]; |
---|
| 28 | newpix[wcs->lat] += pix[1]; |
---|
| 29 | newpix[wcs->spec]+= pix[2]; |
---|
[22] | 30 | |
---|
[128] | 31 | int *stat = new int[npts]; |
---|
[160] | 32 | double *imgcrd = new double[naxis*npts]; |
---|
| 33 | double *tempworld = new double[naxis*npts]; |
---|
[128] | 34 | double *phi = new double[npts]; |
---|
| 35 | double *theta = new double[npts]; |
---|
[160] | 36 | if(flag=wcsp2s(wcs, npts, naxis, newpix, imgcrd, phi, |
---|
| 37 | theta, tempworld, stat)>0){ |
---|
[120] | 38 | std::stringstream errmsg; |
---|
| 39 | errmsg << "WCS Error Code = " <<flag<<": " << wcs_errmsg[flag] |
---|
| 40 | << "\nstat value is " << stat[0] << std::endl; |
---|
| 41 | duchampError("pixToWCSSingle",errmsg.str()); |
---|
[3] | 42 | } |
---|
[128] | 43 | |
---|
| 44 | //return just the spatial/velocity information |
---|
[160] | 45 | world[0] = tempworld[wcs->lng]; |
---|
| 46 | world[1] = tempworld[wcs->lat]; |
---|
| 47 | world[2] = tempworld[wcs->spec]; |
---|
[128] | 48 | |
---|
[3] | 49 | delete [] stat; |
---|
| 50 | delete [] imgcrd; |
---|
[130] | 51 | delete [] tempworld; |
---|
[3] | 52 | delete [] phi; |
---|
| 53 | delete [] theta; |
---|
[103] | 54 | delete [] newpix; |
---|
[3] | 55 | return flag; |
---|
| 56 | } |
---|
| 57 | |
---|
| 58 | int wcsToPixSingle(wcsprm *wcs, const double *world, double *pix) |
---|
| 59 | { |
---|
[22] | 60 | /** |
---|
| 61 | * wcsToPixSingle(wcsprm *wcs, const double *world, double *pix) |
---|
[160] | 62 | * Uses wcs to convert the three-dimensional world coordinate position |
---|
| 63 | * referenced by world to pixel coordinates, which are placed in the |
---|
| 64 | * array pix[]. |
---|
| 65 | * Assumes that world only has one point (three values eg RA Dec Velocity) |
---|
| 66 | * Offsets the pixel values by 1 to account for the C arrays being |
---|
| 67 | * indexed to 0. |
---|
[22] | 68 | */ |
---|
| 69 | |
---|
[160] | 70 | int naxis=wcs->naxis,npts=1,flag; |
---|
[128] | 71 | |
---|
[160] | 72 | double *tempworld = new double[naxis*npts]; |
---|
| 73 | for(int i=0;i<naxis;i++) tempworld[i] = wcs->crval[i]; |
---|
| 74 | tempworld[wcs->lng] = world[0]; |
---|
| 75 | tempworld[wcs->lat] = world[1]; |
---|
| 76 | tempworld[wcs->spec] = world[2]; |
---|
[128] | 77 | |
---|
[160] | 78 | int *stat = new int[npts]; |
---|
| 79 | double *temppix = new double[naxis*npts]; |
---|
| 80 | double *imgcrd = new double[naxis*npts]; |
---|
| 81 | double *phi = new double[npts]; |
---|
| 82 | double *theta = new double[npts]; |
---|
[128] | 83 | |
---|
[160] | 84 | flag=wcss2p(wcs, npts, naxis, tempworld, phi, theta, imgcrd, temppix, stat); |
---|
| 85 | |
---|
| 86 | if( flag > 0 ){ |
---|
[120] | 87 | std::stringstream errmsg; |
---|
| 88 | errmsg << "WCS Error Code = " <<flag<<": " << wcs_errmsg[flag] |
---|
| 89 | << "\nstat value is " << stat[0] << std::endl; |
---|
| 90 | duchampError("wcsToPixSingle",errmsg.str()); |
---|
[3] | 91 | } |
---|
[22] | 92 | |
---|
[160] | 93 | pix[0] = temppix[wcs->lng] - 1.; |
---|
| 94 | pix[1] = temppix[wcs->lat] - 1.; |
---|
| 95 | pix[2] = temppix[wcs->spec] - 1.; |
---|
[128] | 96 | // correct from 1-indexed to 0-indexed pixel array |
---|
| 97 | // and only return the spatial & frequency information |
---|
[22] | 98 | |
---|
[3] | 99 | delete [] stat; |
---|
| 100 | delete [] imgcrd; |
---|
[130] | 101 | delete [] temppix; |
---|
[3] | 102 | delete [] phi; |
---|
| 103 | delete [] theta; |
---|
[130] | 104 | delete [] tempworld; |
---|
[3] | 105 | return flag; |
---|
| 106 | } |
---|
| 107 | |
---|
[22] | 108 | int pixToWCSMulti(wcsprm *wcs, const double *pix, double *world, const int npts) |
---|
| 109 | { |
---|
| 110 | /** |
---|
| 111 | * pixToWCSSingle(wcsprm *wcs, const double *pix, double *world) |
---|
[160] | 112 | * Uses wcs to convert the three-dimensional pixel positions referenced |
---|
| 113 | * by pix to world coordinates, which are placed in the array world[]. |
---|
[22] | 114 | * pix is assumed to hold the positions of npts points. |
---|
[160] | 115 | * Offsets these pixel values by 1 to account for the C arrays being |
---|
| 116 | * indexed to 0. |
---|
[22] | 117 | */ |
---|
| 118 | |
---|
[160] | 119 | int naxis=wcs->naxis,flag; |
---|
[22] | 120 | |
---|
| 121 | // correct from 0-indexed to 1-indexed pixel array |
---|
[128] | 122 | // Add entries for any other axes that are present, keeping the |
---|
| 123 | // order of pixel positions the same |
---|
[160] | 124 | double *newpix = new double[naxis*npts]; |
---|
| 125 | for(int pt=0;pt<npts;pt++){ |
---|
| 126 | for(int i=0;i<naxis;i++) newpix[pt*naxis+i] = 1.; |
---|
| 127 | newpix[pt*naxis+wcs->lng] += pix[pt*3+0]; |
---|
| 128 | newpix[pt*naxis+wcs->lat] += pix[pt*3+1]; |
---|
| 129 | newpix[pt*naxis+wcs->spec] += pix[pt*3+2]; |
---|
[128] | 130 | } |
---|
| 131 | |
---|
| 132 | int *stat = new int[npts]; |
---|
[160] | 133 | double *imgcrd = new double[naxis*npts]; |
---|
| 134 | double *tempworld = new double[naxis*npts]; |
---|
[128] | 135 | double *phi = new double[npts]; |
---|
| 136 | double *theta = new double[npts]; |
---|
| 137 | |
---|
[160] | 138 | if(flag=wcsp2s(wcs, npts, naxis, newpix, imgcrd, |
---|
| 139 | phi, theta, tempworld, stat) >0){ |
---|
[120] | 140 | std::stringstream errmsg; |
---|
| 141 | errmsg << "WCS Error Code = " <<flag<<": " << wcs_errmsg[flag] |
---|
| 142 | << "\nstat value is " << stat[0] << std::endl; |
---|
| 143 | duchampError("pixToWCSMulti",errmsg.str()); |
---|
[22] | 144 | } |
---|
[128] | 145 | else{ |
---|
| 146 | //return just the spatial/velocity information, keeping the |
---|
| 147 | // order of the pixel positions the same. |
---|
[160] | 148 | for(int pt=0;pt<npts;pt++){ |
---|
| 149 | world[pt*3+0] = tempworld[pt*naxis + wcs->lng]; |
---|
| 150 | world[pt*3+1] = tempworld[pt*naxis + wcs->lat]; |
---|
| 151 | world[pt*3+2] = tempworld[pt*naxis + wcs->spec]; |
---|
[128] | 152 | } |
---|
[160] | 153 | |
---|
[128] | 154 | } |
---|
| 155 | |
---|
[22] | 156 | delete [] stat; |
---|
| 157 | delete [] imgcrd; |
---|
[130] | 158 | delete [] tempworld; |
---|
[22] | 159 | delete [] phi; |
---|
| 160 | delete [] theta; |
---|
[38] | 161 | delete [] newpix; |
---|
[22] | 162 | return flag; |
---|
| 163 | } |
---|
| 164 | |
---|
| 165 | int wcsToPixMulti(wcsprm *wcs, const double *world, double *pix, const int npts) |
---|
| 166 | { |
---|
| 167 | /** |
---|
| 168 | * wcsToPixSingle(wcsprm *wcs, const double *world, double *pix) |
---|
[160] | 169 | * Uses wcs to convert the three-dimensional world coordinate position |
---|
| 170 | * referenced by world to pixel coordinates, which are placed in the |
---|
| 171 | * array pix[]. |
---|
[22] | 172 | * world is assumed to hold the positions of npts points. |
---|
[160] | 173 | * Offsets the pixel values by 1 to account for the C arrays being |
---|
| 174 | * indexed to 0. |
---|
[22] | 175 | */ |
---|
[103] | 176 | |
---|
[160] | 177 | int naxis=wcs->naxis,flag=0; |
---|
[128] | 178 | // Test to see if there are other axes present, eg. stokes |
---|
[160] | 179 | if(wcs->naxis>naxis) naxis = wcs->naxis; |
---|
[128] | 180 | |
---|
| 181 | // Add entries for any other axes that are present, keeping the |
---|
| 182 | // order of pixel positions the same |
---|
[160] | 183 | double *tempworld = new double[naxis*npts]; |
---|
| 184 | for(int pt=0;pt<npts;pt++){ |
---|
| 185 | for(int axis=0;axis<naxis;axis++) |
---|
| 186 | tempworld[pt*naxis+axis] = wcs->crval[axis]; |
---|
| 187 | tempworld[pt*naxis + wcs->lng] = world[pt*3+0]; |
---|
| 188 | tempworld[pt*naxis + wcs->lat] = world[pt*3+1]; |
---|
| 189 | tempworld[pt*naxis + wcs->spec] = world[pt*3+2]; |
---|
[128] | 190 | } |
---|
| 191 | |
---|
[22] | 192 | int *stat = new int[npts]; |
---|
[160] | 193 | double *temppix = new double[naxis*npts]; |
---|
| 194 | double *imgcrd = new double[naxis*npts]; |
---|
[22] | 195 | double *phi = new double[npts]; |
---|
| 196 | double *theta = new double[npts]; |
---|
[160] | 197 | if(flag=wcss2p(wcs, npts, naxis, tempworld, |
---|
| 198 | phi, theta, imgcrd, temppix, stat)>0){ |
---|
[120] | 199 | std::stringstream errmsg; |
---|
| 200 | errmsg << "WCS Error Code = " <<flag<<": " <<wcs_errmsg[flag] |
---|
| 201 | << "\nstat value is " << stat[0] << std::endl; |
---|
| 202 | duchampError("wcsToPixMulti",errmsg.str()); |
---|
[22] | 203 | } |
---|
[128] | 204 | else{ |
---|
| 205 | // correct from 1-indexed to 0-indexed pixel array |
---|
| 206 | // and return just the spatial/velocity information, |
---|
| 207 | // keeping the order of the pixel positions the same. |
---|
[160] | 208 | for(int pt=0;pt<npts;pt++){ |
---|
| 209 | pix[pt*naxis + 0] = temppix[pt*naxis + wcs->lng] - 1.; |
---|
| 210 | pix[pt*naxis + 1] = temppix[pt*naxis + wcs->lat] - 1.; |
---|
| 211 | pix[pt*naxis + 2] = temppix[pt*naxis + wcs->spec] - 1.; |
---|
[128] | 212 | } |
---|
| 213 | } |
---|
[22] | 214 | |
---|
| 215 | delete [] stat; |
---|
| 216 | delete [] imgcrd; |
---|
[130] | 217 | delete [] temppix; |
---|
[22] | 218 | delete [] phi; |
---|
| 219 | delete [] theta; |
---|
[130] | 220 | delete [] tempworld; |
---|
[22] | 221 | return flag; |
---|
| 222 | } |
---|
| 223 | |
---|
[160] | 224 | double pixelToVelocity(wcsprm *wcs, double &x, double &y, double &z, |
---|
| 225 | string velUnits) |
---|
[3] | 226 | { |
---|
[22] | 227 | /** |
---|
[160] | 228 | * pixelToVelocity(wcsprm *wcs, double &x, double &y, double &z, |
---|
| 229 | * string velUnits) |
---|
[103] | 230 | * Uses wcs to convert the three-dimensional pixel position (x,y,z) |
---|
| 231 | * to world coordinates. |
---|
| 232 | * Returns the velocity in the units given by velUnits. |
---|
| 233 | */ |
---|
| 234 | |
---|
[160] | 235 | int naxis=wcs->naxis,flag=0; |
---|
[103] | 236 | int *stat = new int[1]; |
---|
[160] | 237 | double *pixcrd = new double[naxis]; |
---|
| 238 | double *imgcrd = new double[naxis]; |
---|
| 239 | double *world = new double[naxis]; |
---|
[103] | 240 | double *phi = new double[1]; |
---|
| 241 | double *theta = new double[1]; |
---|
| 242 | // correct from 0-indexed to 1-indexed pixel array by adding 1 to pixel values. |
---|
[160] | 243 | for(int i=0;i<naxis;i++) pixcrd[i] = 1.; |
---|
| 244 | pixcrd[wcs->lng] += x; |
---|
| 245 | pixcrd[wcs->lat] += y; |
---|
| 246 | pixcrd[wcs->spec]+= z; |
---|
| 247 | |
---|
| 248 | if(flag=wcsp2s(wcs, 1, naxis, pixcrd, imgcrd, phi, theta, world, stat)>0){ |
---|
[120] | 249 | std::stringstream errmsg; |
---|
| 250 | errmsg <<"WCS Error Code = "<<flag<<": "<<wcs_errmsg[flag] |
---|
| 251 | << "\nstat value is "<<stat[0]<<std::endl; |
---|
| 252 | duchampError("pixelToVelocity",errmsg.str()); |
---|
[103] | 253 | } |
---|
| 254 | |
---|
[160] | 255 | double vel = coordToVel(wcs, world[wcs->spec], velUnits); |
---|
[103] | 256 | |
---|
| 257 | delete [] stat; |
---|
| 258 | delete [] pixcrd; |
---|
| 259 | delete [] imgcrd; |
---|
| 260 | delete [] world; |
---|
| 261 | delete [] phi; |
---|
| 262 | delete [] theta; |
---|
| 263 | |
---|
| 264 | return vel; |
---|
| 265 | } |
---|
| 266 | |
---|
| 267 | double coordToVel(wcsprm *wcs, const double coord, string outputUnits) |
---|
| 268 | { |
---|
| 269 | /** |
---|
| 270 | * coordToVel(wcsprm *wcs, const double coord, string outputUnits) |
---|
[22] | 271 | * Convert the wcs coordinate given by coord to a velocity in km/s. |
---|
[160] | 272 | * Does this by checking the ztype parameter in wcs to see if it is |
---|
| 273 | * FREQ or otherwise, and converts accordingly. |
---|
[22] | 274 | */ |
---|
| 275 | |
---|
[120] | 276 | static int errflag = 0; |
---|
[103] | 277 | double scale, offset, power; |
---|
[160] | 278 | int specIndex = wcs->spec; |
---|
| 279 | int flag = wcsunits( wcs->cunit[specIndex], outputUnits.c_str(), |
---|
| 280 | &scale, &offset, &power); |
---|
[103] | 281 | |
---|
| 282 | if(flag>0){ |
---|
[120] | 283 | if(errflag==0){ |
---|
| 284 | std::stringstream errmsg; |
---|
| 285 | errmsg << "WCSUNITS Error Code = " << flag << ":" |
---|
| 286 | << wcsunits_errmsg[flag]; |
---|
[160] | 287 | if(flag==10) errmsg << "\nTried to get conversion from " |
---|
| 288 | << wcs->cunit[specIndex] |
---|
[120] | 289 | << " to " << outputUnits.c_str(); |
---|
| 290 | errmsg << "\nUsing coordinate value instead.\n"; |
---|
| 291 | duchampError("coordToVel", errmsg.str()); |
---|
| 292 | errflag = 1; |
---|
[103] | 293 | } |
---|
| 294 | return coord; |
---|
[3] | 295 | } |
---|
[103] | 296 | else return pow( (coord*scale + offset), power); |
---|
| 297 | |
---|
[3] | 298 | } |
---|
| 299 | |
---|
[103] | 300 | double velToCoord(wcsprm *wcs, const float velocity, string outputUnits) |
---|
[3] | 301 | { |
---|
[22] | 302 | /** |
---|
| 303 | * velToCoord(wcsprm *wcs, const double coord) |
---|
| 304 | * Convert the velocity given to the appropriate world coordinate for wcs. |
---|
[160] | 305 | * Does this by checking the ztype parameter in wcs to see if it is |
---|
| 306 | * FREQ or otherwise, and converts accordingly. |
---|
[22] | 307 | */ |
---|
| 308 | |
---|
[103] | 309 | double scale, offset, power; |
---|
[160] | 310 | int specIndex = wcs->spec; |
---|
| 311 | int flag = wcsunits( wcs->cunit[specIndex], outputUnits.c_str(), |
---|
| 312 | &scale, &offset, &power); |
---|
[103] | 313 | |
---|
| 314 | if(flag>0){ |
---|
[120] | 315 | std::stringstream errmsg; |
---|
| 316 | errmsg << "WCSUNITS Error Code = " << flag << ":" |
---|
| 317 | << wcsunits_errmsg[flag] << "\nUsing coordinate value instead.\n"; |
---|
| 318 | duchampError("velToCoord",errmsg.str()); |
---|
[103] | 319 | return velocity; |
---|
| 320 | } |
---|
| 321 | else return (pow(velocity, 1./power) - offset) / scale; |
---|
[3] | 322 | |
---|
| 323 | } |
---|