[301] | 1 | // ----------------------------------------------------------------------- |
---|
| 2 | // wcsFunctions.cc: Functions to convert between pixel and world |
---|
| 3 | // coordinates -- user-friendly front ends to the |
---|
| 4 | // WCSLIB functions. |
---|
| 5 | // ----------------------------------------------------------------------- |
---|
| 6 | // Copyright (C) 2006, Matthew Whiting, ATNF |
---|
| 7 | // |
---|
| 8 | // This program is free software; you can redistribute it and/or modify it |
---|
| 9 | // under the terms of the GNU General Public License as published by the |
---|
| 10 | // Free Software Foundation; either version 2 of the License, or (at your |
---|
| 11 | // option) any later version. |
---|
| 12 | // |
---|
| 13 | // Duchamp is distributed in the hope that it will be useful, but WITHOUT |
---|
| 14 | // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or |
---|
| 15 | // FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
---|
| 16 | // for more details. |
---|
| 17 | // |
---|
| 18 | // You should have received a copy of the GNU General Public License |
---|
| 19 | // along with Duchamp; if not, write to the Free Software Foundation, |
---|
| 20 | // Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA |
---|
| 21 | // |
---|
| 22 | // Correspondence concerning Duchamp may be directed to: |
---|
| 23 | // Internet email: Matthew.Whiting [at] atnf.csiro.au |
---|
| 24 | // Postal address: Dr. Matthew Whiting |
---|
| 25 | // Australia Telescope National Facility, CSIRO |
---|
| 26 | // PO Box 76 |
---|
| 27 | // Epping NSW 1710 |
---|
| 28 | // AUSTRALIA |
---|
| 29 | // ----------------------------------------------------------------------- |
---|
[3] | 30 | #include <iostream> |
---|
[120] | 31 | #include <sstream> |
---|
[103] | 32 | #include <math.h> |
---|
[394] | 33 | #include <wcslib/wcs.h> |
---|
| 34 | #include <wcslib/wcsunits.h> |
---|
[393] | 35 | #include <duchamp/duchamp.hh> |
---|
| 36 | #include <duchamp/Utils/utils.hh> |
---|
[3] | 37 | |
---|
[204] | 38 | int pixToWCSSingle(struct wcsprm *wcs, const double *pix, double *world) |
---|
[3] | 39 | { |
---|
[22] | 40 | /** |
---|
[160] | 41 | * Uses wcs to convert the three-dimensional pixel position referenced |
---|
| 42 | * by pix to world coordinates, which are placed in the array world[]. |
---|
[22] | 43 | * Assumes that pix only has one point with an x,y,and z pixel positions. |
---|
[160] | 44 | * If there are any other axes present (eg. Stokes), these are set to the |
---|
| 45 | * first pixel in that axis. |
---|
| 46 | * Offsets these pixel values by 1 to account for the C arrays being |
---|
| 47 | * indexed to 0. |
---|
[301] | 48 | * |
---|
| 49 | * \param wcs The wcsprm struct containing the WCS information. |
---|
| 50 | * \param pix The array of pixel coordinates. |
---|
| 51 | * \param world The returned array of world coordinates. |
---|
[22] | 52 | */ |
---|
| 53 | |
---|
[205] | 54 | int naxis=wcs->naxis,npts=1,status; |
---|
[22] | 55 | |
---|
[160] | 56 | double *newpix = new double[naxis*npts]; |
---|
[22] | 57 | // correct from 0-indexed to 1-indexed pixel array |
---|
[160] | 58 | for(int i=0;i<naxis;i++) newpix[i] = 1.; |
---|
| 59 | newpix[wcs->lng] += pix[0]; |
---|
| 60 | newpix[wcs->lat] += pix[1]; |
---|
| 61 | newpix[wcs->spec]+= pix[2]; |
---|
[22] | 62 | |
---|
[128] | 63 | int *stat = new int[npts]; |
---|
[160] | 64 | double *imgcrd = new double[naxis*npts]; |
---|
| 65 | double *tempworld = new double[naxis*npts]; |
---|
[128] | 66 | double *phi = new double[npts]; |
---|
| 67 | double *theta = new double[npts]; |
---|
[301] | 68 | status=wcsp2s(wcs, npts, naxis, newpix, imgcrd, |
---|
| 69 | phi, theta, tempworld, stat); |
---|
[270] | 70 | if(status>0){ |
---|
[120] | 71 | std::stringstream errmsg; |
---|
[205] | 72 | errmsg << "WCS Error Code = " << status <<": " << wcs_errmsg[status] |
---|
[120] | 73 | << "\nstat value is " << stat[0] << std::endl; |
---|
[378] | 74 | duchamp::duchampError("pixToWCSSingle",errmsg.str()); |
---|
[3] | 75 | } |
---|
[128] | 76 | |
---|
| 77 | //return just the spatial/velocity information |
---|
[160] | 78 | world[0] = tempworld[wcs->lng]; |
---|
| 79 | world[1] = tempworld[wcs->lat]; |
---|
| 80 | world[2] = tempworld[wcs->spec]; |
---|
[128] | 81 | |
---|
[3] | 82 | delete [] stat; |
---|
| 83 | delete [] imgcrd; |
---|
[130] | 84 | delete [] tempworld; |
---|
[3] | 85 | delete [] phi; |
---|
| 86 | delete [] theta; |
---|
[103] | 87 | delete [] newpix; |
---|
[205] | 88 | return status; |
---|
[3] | 89 | } |
---|
| 90 | |
---|
[204] | 91 | int wcsToPixSingle(struct wcsprm *wcs, const double *world, double *pix) |
---|
[3] | 92 | { |
---|
[22] | 93 | /** |
---|
[160] | 94 | * Uses wcs to convert the three-dimensional world coordinate position |
---|
| 95 | * referenced by world to pixel coordinates, which are placed in the |
---|
| 96 | * array pix[]. |
---|
| 97 | * Assumes that world only has one point (three values eg RA Dec Velocity) |
---|
| 98 | * Offsets the pixel values by 1 to account for the C arrays being |
---|
| 99 | * indexed to 0. |
---|
[301] | 100 | * |
---|
| 101 | * \param wcs The wcsprm struct containing the WCS information. |
---|
| 102 | * \param world The array of world coordinates. |
---|
| 103 | * \param pix The returned array of pixel coordinates. |
---|
[22] | 104 | */ |
---|
| 105 | |
---|
[205] | 106 | int naxis=wcs->naxis,npts=1,status; |
---|
[128] | 107 | |
---|
[160] | 108 | double *tempworld = new double[naxis*npts]; |
---|
| 109 | for(int i=0;i<naxis;i++) tempworld[i] = wcs->crval[i]; |
---|
| 110 | tempworld[wcs->lng] = world[0]; |
---|
| 111 | tempworld[wcs->lat] = world[1]; |
---|
| 112 | tempworld[wcs->spec] = world[2]; |
---|
[128] | 113 | |
---|
[160] | 114 | int *stat = new int[npts]; |
---|
| 115 | double *temppix = new double[naxis*npts]; |
---|
| 116 | double *imgcrd = new double[naxis*npts]; |
---|
| 117 | double *phi = new double[npts]; |
---|
| 118 | double *theta = new double[npts]; |
---|
[128] | 119 | |
---|
[301] | 120 | status=wcss2p(wcs, npts, naxis, tempworld, |
---|
| 121 | phi, theta, imgcrd, temppix, stat); |
---|
[160] | 122 | |
---|
[205] | 123 | if( status > 0 ){ |
---|
[120] | 124 | std::stringstream errmsg; |
---|
[205] | 125 | errmsg << "WCS Error Code = " <<status<<": " << wcs_errmsg[status] |
---|
[120] | 126 | << "\nstat value is " << stat[0] << std::endl; |
---|
[378] | 127 | duchamp::duchampError("wcsToPixSingle",errmsg.str()); |
---|
[3] | 128 | } |
---|
[22] | 129 | |
---|
[160] | 130 | pix[0] = temppix[wcs->lng] - 1.; |
---|
| 131 | pix[1] = temppix[wcs->lat] - 1.; |
---|
| 132 | pix[2] = temppix[wcs->spec] - 1.; |
---|
[128] | 133 | // correct from 1-indexed to 0-indexed pixel array |
---|
| 134 | // and only return the spatial & frequency information |
---|
[22] | 135 | |
---|
[3] | 136 | delete [] stat; |
---|
| 137 | delete [] imgcrd; |
---|
[130] | 138 | delete [] temppix; |
---|
[3] | 139 | delete [] phi; |
---|
| 140 | delete [] theta; |
---|
[130] | 141 | delete [] tempworld; |
---|
[205] | 142 | return status; |
---|
[3] | 143 | } |
---|
| 144 | |
---|
[301] | 145 | int pixToWCSMulti(struct wcsprm *wcs, const double *pix, |
---|
| 146 | double *world, const int npts) |
---|
[22] | 147 | { |
---|
| 148 | /** |
---|
[160] | 149 | * Uses wcs to convert the three-dimensional pixel positions referenced |
---|
| 150 | * by pix to world coordinates, which are placed in the array world[]. |
---|
[22] | 151 | * pix is assumed to hold the positions of npts points. |
---|
[160] | 152 | * Offsets these pixel values by 1 to account for the C arrays being |
---|
| 153 | * indexed to 0. |
---|
[301] | 154 | * |
---|
| 155 | * \param wcs The wcsprm struct containing the WCS information. |
---|
| 156 | * \param pix The array of pixel coordinates. |
---|
| 157 | * \param world The returned array of world coordinates. |
---|
| 158 | * \param npts The number of distinct pixels in the arrays. |
---|
[22] | 159 | */ |
---|
| 160 | |
---|
[205] | 161 | int naxis=wcs->naxis,status; |
---|
[22] | 162 | |
---|
| 163 | // correct from 0-indexed to 1-indexed pixel array |
---|
[128] | 164 | // Add entries for any other axes that are present, keeping the |
---|
| 165 | // order of pixel positions the same |
---|
[160] | 166 | double *newpix = new double[naxis*npts]; |
---|
| 167 | for(int pt=0;pt<npts;pt++){ |
---|
| 168 | for(int i=0;i<naxis;i++) newpix[pt*naxis+i] = 1.; |
---|
| 169 | newpix[pt*naxis+wcs->lng] += pix[pt*3+0]; |
---|
| 170 | newpix[pt*naxis+wcs->lat] += pix[pt*3+1]; |
---|
| 171 | newpix[pt*naxis+wcs->spec] += pix[pt*3+2]; |
---|
[128] | 172 | } |
---|
| 173 | |
---|
| 174 | int *stat = new int[npts]; |
---|
[160] | 175 | double *imgcrd = new double[naxis*npts]; |
---|
| 176 | double *tempworld = new double[naxis*npts]; |
---|
[128] | 177 | double *phi = new double[npts]; |
---|
| 178 | double *theta = new double[npts]; |
---|
[301] | 179 | status=wcsp2s(wcs, npts, naxis, newpix, imgcrd, |
---|
| 180 | phi, theta, tempworld, stat); |
---|
[270] | 181 | if(status>0){ |
---|
[120] | 182 | std::stringstream errmsg; |
---|
[205] | 183 | errmsg << "WCS Error Code = " <<status<<": " << wcs_errmsg[status] |
---|
[120] | 184 | << "\nstat value is " << stat[0] << std::endl; |
---|
[378] | 185 | duchamp::duchampError("pixToWCSMulti",errmsg.str()); |
---|
[22] | 186 | } |
---|
[128] | 187 | else{ |
---|
| 188 | //return just the spatial/velocity information, keeping the |
---|
| 189 | // order of the pixel positions the same. |
---|
[160] | 190 | for(int pt=0;pt<npts;pt++){ |
---|
| 191 | world[pt*3+0] = tempworld[pt*naxis + wcs->lng]; |
---|
| 192 | world[pt*3+1] = tempworld[pt*naxis + wcs->lat]; |
---|
| 193 | world[pt*3+2] = tempworld[pt*naxis + wcs->spec]; |
---|
[128] | 194 | } |
---|
[160] | 195 | |
---|
[128] | 196 | } |
---|
| 197 | |
---|
[22] | 198 | delete [] stat; |
---|
| 199 | delete [] imgcrd; |
---|
[130] | 200 | delete [] tempworld; |
---|
[22] | 201 | delete [] phi; |
---|
| 202 | delete [] theta; |
---|
[38] | 203 | delete [] newpix; |
---|
[205] | 204 | return status; |
---|
[22] | 205 | } |
---|
| 206 | |
---|
[301] | 207 | int wcsToPixMulti(struct wcsprm *wcs, const double *world, |
---|
| 208 | double *pix, const int npts) |
---|
[22] | 209 | { |
---|
| 210 | /** |
---|
[160] | 211 | * Uses wcs to convert the three-dimensional world coordinate position |
---|
| 212 | * referenced by world to pixel coordinates, which are placed in the |
---|
| 213 | * array pix[]. |
---|
[22] | 214 | * world is assumed to hold the positions of npts points. |
---|
[160] | 215 | * Offsets the pixel values by 1 to account for the C arrays being |
---|
| 216 | * indexed to 0. |
---|
[301] | 217 | * |
---|
| 218 | * \param wcs The wcsprm struct containing the WCS information. |
---|
| 219 | * \param world The array of world coordinates. |
---|
| 220 | * \param pix The returned array of pixel coordinates. |
---|
| 221 | * \param npts The number of distinct pixels in the arrays. |
---|
[22] | 222 | */ |
---|
[103] | 223 | |
---|
[205] | 224 | int naxis=wcs->naxis,status=0; |
---|
[128] | 225 | // Test to see if there are other axes present, eg. stokes |
---|
[160] | 226 | if(wcs->naxis>naxis) naxis = wcs->naxis; |
---|
[128] | 227 | |
---|
| 228 | // Add entries for any other axes that are present, keeping the |
---|
| 229 | // order of pixel positions the same |
---|
[160] | 230 | double *tempworld = new double[naxis*npts]; |
---|
| 231 | for(int pt=0;pt<npts;pt++){ |
---|
| 232 | for(int axis=0;axis<naxis;axis++) |
---|
| 233 | tempworld[pt*naxis+axis] = wcs->crval[axis]; |
---|
| 234 | tempworld[pt*naxis + wcs->lng] = world[pt*3+0]; |
---|
| 235 | tempworld[pt*naxis + wcs->lat] = world[pt*3+1]; |
---|
| 236 | tempworld[pt*naxis + wcs->spec] = world[pt*3+2]; |
---|
[128] | 237 | } |
---|
| 238 | |
---|
[22] | 239 | int *stat = new int[npts]; |
---|
[160] | 240 | double *temppix = new double[naxis*npts]; |
---|
| 241 | double *imgcrd = new double[naxis*npts]; |
---|
[22] | 242 | double *phi = new double[npts]; |
---|
| 243 | double *theta = new double[npts]; |
---|
[301] | 244 | status=wcss2p(wcs,npts,naxis,tempworld,phi,theta, |
---|
| 245 | imgcrd,temppix,stat); |
---|
[270] | 246 | if(status>0){ |
---|
[120] | 247 | std::stringstream errmsg; |
---|
[205] | 248 | errmsg << "WCS Error Code = " <<status<<": " <<wcs_errmsg[status] |
---|
[120] | 249 | << "\nstat value is " << stat[0] << std::endl; |
---|
[378] | 250 | duchamp::duchampError("wcsToPixMulti",errmsg.str()); |
---|
[22] | 251 | } |
---|
[128] | 252 | else{ |
---|
| 253 | // correct from 1-indexed to 0-indexed pixel array |
---|
| 254 | // and return just the spatial/velocity information, |
---|
| 255 | // keeping the order of the pixel positions the same. |
---|
[160] | 256 | for(int pt=0;pt<npts;pt++){ |
---|
| 257 | pix[pt*naxis + 0] = temppix[pt*naxis + wcs->lng] - 1.; |
---|
| 258 | pix[pt*naxis + 1] = temppix[pt*naxis + wcs->lat] - 1.; |
---|
| 259 | pix[pt*naxis + 2] = temppix[pt*naxis + wcs->spec] - 1.; |
---|
[128] | 260 | } |
---|
| 261 | } |
---|
[22] | 262 | |
---|
| 263 | delete [] stat; |
---|
| 264 | delete [] imgcrd; |
---|
[130] | 265 | delete [] temppix; |
---|
[22] | 266 | delete [] phi; |
---|
| 267 | delete [] theta; |
---|
[130] | 268 | delete [] tempworld; |
---|
[205] | 269 | return status; |
---|
[22] | 270 | } |
---|
| 271 | |
---|
[204] | 272 | double pixelToVelocity(struct wcsprm *wcs, double &x, double &y, double &z, |
---|
[232] | 273 | std::string velUnits) |
---|
[3] | 274 | { |
---|
[22] | 275 | /** |
---|
[103] | 276 | * Uses wcs to convert the three-dimensional pixel position (x,y,z) |
---|
| 277 | * to world coordinates. |
---|
| 278 | * Returns the velocity in the units given by velUnits. |
---|
[301] | 279 | * |
---|
| 280 | * \param wcs The wcsprm struct containing the WCS information. |
---|
| 281 | * \param x The x pixel coordinate. |
---|
| 282 | * \param y The y pixel coordinate. |
---|
| 283 | * \param z The z pixel coordinate. |
---|
| 284 | * \param velUnits The string containing the units of the output velocity. |
---|
| 285 | * \return The velocity of the pixel. |
---|
[103] | 286 | */ |
---|
| 287 | |
---|
[205] | 288 | int naxis=wcs->naxis,status=0; |
---|
[103] | 289 | int *stat = new int[1]; |
---|
[160] | 290 | double *pixcrd = new double[naxis]; |
---|
| 291 | double *imgcrd = new double[naxis]; |
---|
| 292 | double *world = new double[naxis]; |
---|
[103] | 293 | double *phi = new double[1]; |
---|
| 294 | double *theta = new double[1]; |
---|
[293] | 295 | // correct from 0-indexed to 1-indexed pixel array by adding 1 to |
---|
| 296 | // pixel values. |
---|
[160] | 297 | for(int i=0;i<naxis;i++) pixcrd[i] = 1.; |
---|
| 298 | pixcrd[wcs->lng] += x; |
---|
| 299 | pixcrd[wcs->lat] += y; |
---|
| 300 | pixcrd[wcs->spec]+= z; |
---|
[301] | 301 | status=wcsp2s(wcs, 1, naxis, pixcrd, imgcrd, |
---|
| 302 | phi, theta, world, stat); |
---|
[270] | 303 | if(status>0){ |
---|
[120] | 304 | std::stringstream errmsg; |
---|
[205] | 305 | errmsg <<"WCS Error Code = "<<status<<": "<<wcs_errmsg[status] |
---|
[120] | 306 | << "\nstat value is "<<stat[0]<<std::endl; |
---|
[378] | 307 | duchamp::duchampError("pixelToVelocity",errmsg.str()); |
---|
[103] | 308 | } |
---|
| 309 | |
---|
[160] | 310 | double vel = coordToVel(wcs, world[wcs->spec], velUnits); |
---|
[103] | 311 | |
---|
| 312 | delete [] stat; |
---|
| 313 | delete [] pixcrd; |
---|
| 314 | delete [] imgcrd; |
---|
| 315 | delete [] world; |
---|
| 316 | delete [] phi; |
---|
| 317 | delete [] theta; |
---|
| 318 | |
---|
| 319 | return vel; |
---|
| 320 | } |
---|
| 321 | |
---|
[301] | 322 | double coordToVel(struct wcsprm *wcs, const double coord, |
---|
| 323 | std::string outputUnits) |
---|
[103] | 324 | { |
---|
| 325 | /** |
---|
[22] | 326 | * Convert the wcs coordinate given by coord to a velocity in km/s. |
---|
[160] | 327 | * Does this by checking the ztype parameter in wcs to see if it is |
---|
| 328 | * FREQ or otherwise, and converts accordingly. |
---|
[301] | 329 | * |
---|
| 330 | * \param wcs The wcsprm struct containing the WCS information. |
---|
| 331 | * \param coord The input WCS coordinate. |
---|
| 332 | * \param outputUnits The string containing the units of the output velocity. |
---|
| 333 | * \return The velocity of the input coord value. |
---|
[22] | 334 | */ |
---|
| 335 | |
---|
[120] | 336 | static int errflag = 0; |
---|
[103] | 337 | double scale, offset, power; |
---|
[160] | 338 | int specIndex = wcs->spec; |
---|
[205] | 339 | int status = wcsunits( wcs->cunit[specIndex], outputUnits.c_str(), |
---|
| 340 | &scale, &offset, &power); |
---|
[103] | 341 | |
---|
[205] | 342 | if(status > 0){ |
---|
[120] | 343 | if(errflag==0){ |
---|
| 344 | std::stringstream errmsg; |
---|
[205] | 345 | errmsg << "WCSUNITS Error Code = " << status << ":" |
---|
| 346 | << wcsunits_errmsg[status]; |
---|
| 347 | if(status == 10) errmsg << "\nTried to get conversion from " |
---|
| 348 | << wcs->cunit[specIndex] |
---|
| 349 | << " to " << outputUnits.c_str(); |
---|
[120] | 350 | errmsg << "\nUsing coordinate value instead.\n"; |
---|
[378] | 351 | duchamp::duchampError("coordToVel", errmsg.str()); |
---|
[120] | 352 | errflag = 1; |
---|
[103] | 353 | } |
---|
| 354 | return coord; |
---|
[3] | 355 | } |
---|
[103] | 356 | else return pow( (coord*scale + offset), power); |
---|
| 357 | |
---|
[3] | 358 | } |
---|
| 359 | |
---|
[301] | 360 | double velToCoord(struct wcsprm *wcs, const float velocity, |
---|
| 361 | std::string outputUnits) |
---|
[3] | 362 | { |
---|
[22] | 363 | /** |
---|
| 364 | * Convert the velocity given to the appropriate world coordinate for wcs. |
---|
[160] | 365 | * Does this by checking the ztype parameter in wcs to see if it is |
---|
| 366 | * FREQ or otherwise, and converts accordingly. |
---|
[301] | 367 | * |
---|
| 368 | * \param wcs The wcsprm struct containing the WCS information. |
---|
| 369 | * \param velocity The input velocity. |
---|
| 370 | * \param outputUnits The string containing the units of the input velocity. |
---|
| 371 | * \return The WCS coordinate corresponding to the input velocity. |
---|
[22] | 372 | */ |
---|
| 373 | |
---|
[103] | 374 | double scale, offset, power; |
---|
[160] | 375 | int specIndex = wcs->spec; |
---|
[301] | 376 | int status = wcsunits( outputUnits.c_str(), wcs->cunit[specIndex], |
---|
[205] | 377 | &scale, &offset, &power); |
---|
[103] | 378 | |
---|
[205] | 379 | if(status > 0){ |
---|
[120] | 380 | std::stringstream errmsg; |
---|
[205] | 381 | errmsg << "WCSUNITS Error Code = " << status << ":" |
---|
| 382 | << wcsunits_errmsg[status] << "\nUsing coordinate value instead.\n"; |
---|
[378] | 383 | duchamp::duchampError("velToCoord",errmsg.str()); |
---|
[103] | 384 | return velocity; |
---|
| 385 | } |
---|
| 386 | else return (pow(velocity, 1./power) - offset) / scale; |
---|
[3] | 387 | |
---|
| 388 | } |
---|