[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 | |
---|
[913] | 38 | using namespace duchamp; |
---|
| 39 | |
---|
[204] | 40 | int pixToWCSSingle(struct wcsprm *wcs, const double *pix, double *world) |
---|
[3] | 41 | { |
---|
[528] | 42 | /// @details |
---|
| 43 | /// Uses wcs to convert the three-dimensional pixel position referenced |
---|
| 44 | /// by pix to world coordinates, which are placed in the array world[]. |
---|
| 45 | /// Assumes that pix only has one point with an x,y,and z pixel positions. |
---|
| 46 | /// If there are any other axes present (eg. Stokes), these are set to the |
---|
| 47 | /// first pixel in that axis. |
---|
| 48 | /// Offsets these pixel values by 1 to account for the C arrays being |
---|
| 49 | /// indexed to 0. |
---|
| 50 | /// |
---|
| 51 | /// \param wcs The wcsprm struct containing the WCS information. |
---|
| 52 | /// \param pix The array of pixel coordinates. |
---|
| 53 | /// \param world The returned array of world coordinates. |
---|
[22] | 54 | |
---|
[205] | 55 | int naxis=wcs->naxis,npts=1,status; |
---|
[478] | 56 | int specAxis = wcs->spec; |
---|
| 57 | if(specAxis<0) specAxis=2; |
---|
[491] | 58 | if(specAxis>=naxis) specAxis = naxis-1; |
---|
[22] | 59 | |
---|
[160] | 60 | double *newpix = new double[naxis*npts]; |
---|
[22] | 61 | // correct from 0-indexed to 1-indexed pixel array |
---|
[160] | 62 | for(int i=0;i<naxis;i++) newpix[i] = 1.; |
---|
| 63 | newpix[wcs->lng] += pix[0]; |
---|
| 64 | newpix[wcs->lat] += pix[1]; |
---|
[478] | 65 | newpix[specAxis]+= pix[2]; |
---|
[22] | 66 | |
---|
[128] | 67 | int *stat = new int[npts]; |
---|
[160] | 68 | double *imgcrd = new double[naxis*npts]; |
---|
| 69 | double *tempworld = new double[naxis*npts]; |
---|
[128] | 70 | double *phi = new double[npts]; |
---|
| 71 | double *theta = new double[npts]; |
---|
[301] | 72 | status=wcsp2s(wcs, npts, naxis, newpix, imgcrd, |
---|
| 73 | phi, theta, tempworld, stat); |
---|
[270] | 74 | if(status>0){ |
---|
[933] | 75 | DUCHAMPTHROW("pixToWCSSingle","WCS Error Code = " << status <<": stat="<<stat[0] << " : " << wcs_errmsg[status]); |
---|
[3] | 76 | } |
---|
[128] | 77 | |
---|
| 78 | //return just the spatial/velocity information |
---|
[160] | 79 | world[0] = tempworld[wcs->lng]; |
---|
| 80 | world[1] = tempworld[wcs->lat]; |
---|
[478] | 81 | world[2] = tempworld[specAxis]; |
---|
[128] | 82 | |
---|
[3] | 83 | delete [] stat; |
---|
| 84 | delete [] imgcrd; |
---|
[130] | 85 | delete [] tempworld; |
---|
[3] | 86 | delete [] phi; |
---|
| 87 | delete [] theta; |
---|
[103] | 88 | delete [] newpix; |
---|
[205] | 89 | return status; |
---|
[3] | 90 | } |
---|
| 91 | |
---|
[204] | 92 | int wcsToPixSingle(struct wcsprm *wcs, const double *world, double *pix) |
---|
[3] | 93 | { |
---|
[528] | 94 | /// @details |
---|
| 95 | /// Uses wcs to convert the three-dimensional world coordinate position |
---|
| 96 | /// referenced by world to pixel coordinates, which are placed in the |
---|
| 97 | /// array pix[]. |
---|
| 98 | /// Assumes that world only has one point (three values eg RA Dec Velocity) |
---|
| 99 | /// Offsets the pixel values by 1 to account for the C arrays being |
---|
| 100 | /// indexed to 0. |
---|
| 101 | /// |
---|
| 102 | /// \param wcs The wcsprm struct containing the WCS information. |
---|
| 103 | /// \param world The array of world coordinates. |
---|
| 104 | /// \param pix The returned array of pixel coordinates. |
---|
[22] | 105 | |
---|
[205] | 106 | int naxis=wcs->naxis,npts=1,status; |
---|
[478] | 107 | int specAxis = wcs->spec; |
---|
| 108 | if(specAxis<0) specAxis=2; |
---|
[491] | 109 | if(specAxis>=naxis) specAxis = naxis-1; |
---|
[128] | 110 | |
---|
[160] | 111 | double *tempworld = new double[naxis*npts]; |
---|
| 112 | for(int i=0;i<naxis;i++) tempworld[i] = wcs->crval[i]; |
---|
| 113 | tempworld[wcs->lng] = world[0]; |
---|
| 114 | tempworld[wcs->lat] = world[1]; |
---|
[478] | 115 | tempworld[specAxis] = world[2]; |
---|
[128] | 116 | |
---|
[160] | 117 | int *stat = new int[npts]; |
---|
| 118 | double *temppix = new double[naxis*npts]; |
---|
| 119 | double *imgcrd = new double[naxis*npts]; |
---|
| 120 | double *phi = new double[npts]; |
---|
| 121 | double *theta = new double[npts]; |
---|
[128] | 122 | |
---|
[301] | 123 | status=wcss2p(wcs, npts, naxis, tempworld, |
---|
| 124 | phi, theta, imgcrd, temppix, stat); |
---|
[160] | 125 | |
---|
[205] | 126 | if( status > 0 ){ |
---|
[933] | 127 | DUCHAMPTHROW("wcsToPixSingle","WCS Error Code = " << status <<": stat="<<stat[0] << " : " << wcs_errmsg[status]); |
---|
[3] | 128 | } |
---|
[22] | 129 | |
---|
[160] | 130 | pix[0] = temppix[wcs->lng] - 1.; |
---|
| 131 | pix[1] = temppix[wcs->lat] - 1.; |
---|
[478] | 132 | pix[2] = temppix[specAxis] - 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 | { |
---|
[528] | 148 | /// @details |
---|
| 149 | /// Uses wcs to convert the three-dimensional pixel positions referenced |
---|
| 150 | /// by pix to world coordinates, which are placed in the array world[]. |
---|
| 151 | /// pix is assumed to hold the positions of npts points. |
---|
| 152 | /// Offsets these pixel values by 1 to account for the C arrays being |
---|
| 153 | /// indexed to 0. |
---|
| 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 | |
---|
[205] | 160 | int naxis=wcs->naxis,status; |
---|
[478] | 161 | int specAxis = wcs->spec; |
---|
| 162 | if(specAxis<0) specAxis=2; |
---|
[491] | 163 | if(specAxis>=naxis) specAxis = naxis-1; |
---|
[22] | 164 | |
---|
| 165 | // correct from 0-indexed to 1-indexed pixel array |
---|
[128] | 166 | // Add entries for any other axes that are present, keeping the |
---|
| 167 | // order of pixel positions the same |
---|
[160] | 168 | double *newpix = new double[naxis*npts]; |
---|
| 169 | for(int pt=0;pt<npts;pt++){ |
---|
| 170 | for(int i=0;i<naxis;i++) newpix[pt*naxis+i] = 1.; |
---|
| 171 | newpix[pt*naxis+wcs->lng] += pix[pt*3+0]; |
---|
| 172 | newpix[pt*naxis+wcs->lat] += pix[pt*3+1]; |
---|
[478] | 173 | newpix[pt*naxis+specAxis] += pix[pt*3+2]; |
---|
[128] | 174 | } |
---|
| 175 | |
---|
| 176 | int *stat = new int[npts]; |
---|
[160] | 177 | double *imgcrd = new double[naxis*npts]; |
---|
| 178 | double *tempworld = new double[naxis*npts]; |
---|
[128] | 179 | double *phi = new double[npts]; |
---|
| 180 | double *theta = new double[npts]; |
---|
[301] | 181 | status=wcsp2s(wcs, npts, naxis, newpix, imgcrd, |
---|
| 182 | phi, theta, tempworld, stat); |
---|
[270] | 183 | if(status>0){ |
---|
[924] | 184 | std::stringstream statstr; |
---|
| 185 | statstr<<"|"; |
---|
[933] | 186 | for(int i=0;i<npts;i++) statstr<<stat[i]<<"|"; |
---|
| 187 | DUCHAMPTHROW("pixToWCSMulti","WCS Error Code = " << status <<": stat="<<statstr.str() << " : " << wcs_errmsg[status]); |
---|
[22] | 188 | } |
---|
[128] | 189 | else{ |
---|
| 190 | //return just the spatial/velocity information, keeping the |
---|
| 191 | // order of the pixel positions the same. |
---|
[160] | 192 | for(int pt=0;pt<npts;pt++){ |
---|
| 193 | world[pt*3+0] = tempworld[pt*naxis + wcs->lng]; |
---|
| 194 | world[pt*3+1] = tempworld[pt*naxis + wcs->lat]; |
---|
[478] | 195 | world[pt*3+2] = tempworld[pt*naxis + specAxis]; |
---|
[128] | 196 | } |
---|
[160] | 197 | |
---|
[128] | 198 | } |
---|
| 199 | |
---|
[22] | 200 | delete [] stat; |
---|
| 201 | delete [] imgcrd; |
---|
[130] | 202 | delete [] tempworld; |
---|
[22] | 203 | delete [] phi; |
---|
| 204 | delete [] theta; |
---|
[38] | 205 | delete [] newpix; |
---|
[205] | 206 | return status; |
---|
[22] | 207 | } |
---|
| 208 | |
---|
[301] | 209 | int wcsToPixMulti(struct wcsprm *wcs, const double *world, |
---|
| 210 | double *pix, const int npts) |
---|
[22] | 211 | { |
---|
[528] | 212 | /// @details |
---|
| 213 | /// Uses wcs to convert the three-dimensional world coordinate position |
---|
| 214 | /// referenced by world to pixel coordinates, which are placed in the |
---|
| 215 | /// array pix[]. |
---|
| 216 | /// world is assumed to hold the positions of npts points. |
---|
| 217 | /// Offsets the pixel values by 1 to account for the C arrays being |
---|
| 218 | /// indexed to 0. |
---|
| 219 | /// |
---|
| 220 | /// \param wcs The wcsprm struct containing the WCS information. |
---|
| 221 | /// \param world The array of world coordinates. |
---|
| 222 | /// \param pix The returned array of pixel coordinates. |
---|
| 223 | /// \param npts The number of distinct pixels in the arrays. |
---|
[103] | 224 | |
---|
[205] | 225 | int naxis=wcs->naxis,status=0; |
---|
[478] | 226 | int specAxis = wcs->spec; |
---|
| 227 | if(specAxis<0) specAxis=2; |
---|
[491] | 228 | if(specAxis>=naxis) specAxis = naxis-1; |
---|
[478] | 229 | |
---|
[128] | 230 | // Test to see if there are other axes present, eg. stokes |
---|
[160] | 231 | if(wcs->naxis>naxis) naxis = wcs->naxis; |
---|
[128] | 232 | |
---|
| 233 | // Add entries for any other axes that are present, keeping the |
---|
| 234 | // order of pixel positions the same |
---|
[160] | 235 | double *tempworld = new double[naxis*npts]; |
---|
| 236 | for(int pt=0;pt<npts;pt++){ |
---|
| 237 | for(int axis=0;axis<naxis;axis++) |
---|
| 238 | tempworld[pt*naxis+axis] = wcs->crval[axis]; |
---|
| 239 | tempworld[pt*naxis + wcs->lng] = world[pt*3+0]; |
---|
| 240 | tempworld[pt*naxis + wcs->lat] = world[pt*3+1]; |
---|
[478] | 241 | tempworld[pt*naxis + specAxis] = world[pt*3+2]; |
---|
[128] | 242 | } |
---|
| 243 | |
---|
[22] | 244 | int *stat = new int[npts]; |
---|
[160] | 245 | double *temppix = new double[naxis*npts]; |
---|
| 246 | double *imgcrd = new double[naxis*npts]; |
---|
[22] | 247 | double *phi = new double[npts]; |
---|
| 248 | double *theta = new double[npts]; |
---|
[301] | 249 | status=wcss2p(wcs,npts,naxis,tempworld,phi,theta, |
---|
| 250 | imgcrd,temppix,stat); |
---|
[270] | 251 | if(status>0){ |
---|
[924] | 252 | std::stringstream statstr; |
---|
| 253 | statstr<<"|"; |
---|
[933] | 254 | for(int i=0;i<npts;i++) statstr<<stat[i]<<"|"; |
---|
| 255 | DUCHAMPTHROW("wcsToPixMulti","WCS Error Code = " << status <<": stat="<<statstr.str() << " : " << wcs_errmsg[status]); |
---|
[22] | 256 | } |
---|
[128] | 257 | else{ |
---|
| 258 | // correct from 1-indexed to 0-indexed pixel array |
---|
| 259 | // and return just the spatial/velocity information, |
---|
| 260 | // keeping the order of the pixel positions the same. |
---|
[160] | 261 | for(int pt=0;pt<npts;pt++){ |
---|
| 262 | pix[pt*naxis + 0] = temppix[pt*naxis + wcs->lng] - 1.; |
---|
| 263 | pix[pt*naxis + 1] = temppix[pt*naxis + wcs->lat] - 1.; |
---|
[478] | 264 | pix[pt*naxis + 2] = temppix[pt*naxis + specAxis] - 1.; |
---|
[128] | 265 | } |
---|
| 266 | } |
---|
[22] | 267 | |
---|
| 268 | delete [] stat; |
---|
| 269 | delete [] imgcrd; |
---|
[130] | 270 | delete [] temppix; |
---|
[22] | 271 | delete [] phi; |
---|
| 272 | delete [] theta; |
---|
[130] | 273 | delete [] tempworld; |
---|
[205] | 274 | return status; |
---|
[22] | 275 | } |
---|
| 276 | |
---|
[204] | 277 | double pixelToVelocity(struct wcsprm *wcs, double &x, double &y, double &z, |
---|
[232] | 278 | std::string velUnits) |
---|
[3] | 279 | { |
---|
[528] | 280 | /// @details |
---|
| 281 | /// Uses wcs to convert the three-dimensional pixel position (x,y,z) |
---|
| 282 | /// to world coordinates. |
---|
| 283 | /// Returns the velocity in the units given by velUnits. |
---|
| 284 | /// |
---|
| 285 | /// \param wcs The wcsprm struct containing the WCS information. |
---|
| 286 | /// \param x The x pixel coordinate. |
---|
| 287 | /// \param y The y pixel coordinate. |
---|
| 288 | /// \param z The z pixel coordinate. |
---|
| 289 | /// \param velUnits The string containing the units of the output velocity. |
---|
| 290 | /// \return The velocity of the pixel. |
---|
[103] | 291 | |
---|
[205] | 292 | int naxis=wcs->naxis,status=0; |
---|
[103] | 293 | int *stat = new int[1]; |
---|
[160] | 294 | double *pixcrd = new double[naxis]; |
---|
| 295 | double *imgcrd = new double[naxis]; |
---|
| 296 | double *world = new double[naxis]; |
---|
[103] | 297 | double *phi = new double[1]; |
---|
| 298 | double *theta = new double[1]; |
---|
[478] | 299 | |
---|
| 300 | int specAxis = wcs->spec; |
---|
| 301 | if(specAxis<0) specAxis=2; |
---|
[491] | 302 | if(specAxis>=naxis) specAxis = naxis-1; |
---|
[478] | 303 | |
---|
[293] | 304 | // correct from 0-indexed to 1-indexed pixel array by adding 1 to |
---|
| 305 | // pixel values. |
---|
[160] | 306 | for(int i=0;i<naxis;i++) pixcrd[i] = 1.; |
---|
| 307 | pixcrd[wcs->lng] += x; |
---|
| 308 | pixcrd[wcs->lat] += y; |
---|
[478] | 309 | pixcrd[specAxis]+= z; |
---|
[301] | 310 | status=wcsp2s(wcs, 1, naxis, pixcrd, imgcrd, |
---|
| 311 | phi, theta, world, stat); |
---|
[270] | 312 | if(status>0){ |
---|
[924] | 313 | std::stringstream statstr; |
---|
| 314 | statstr<<"|"; |
---|
| 315 | for(int i=0;i<naxis;i++) statstr<<stat[i]<<"|"; |
---|
[933] | 316 | DUCHAMPTHROW("pixelToVelocity","WCS Error Code = " << status <<": stat="<<stat[0] << " : " << wcs_errmsg[status]); |
---|
[103] | 317 | } |
---|
| 318 | |
---|
[478] | 319 | double vel = coordToVel(wcs, world[specAxis], velUnits); |
---|
[103] | 320 | |
---|
| 321 | delete [] stat; |
---|
| 322 | delete [] pixcrd; |
---|
| 323 | delete [] imgcrd; |
---|
| 324 | delete [] world; |
---|
| 325 | delete [] phi; |
---|
| 326 | delete [] theta; |
---|
| 327 | |
---|
| 328 | return vel; |
---|
| 329 | } |
---|
| 330 | |
---|
[301] | 331 | double coordToVel(struct wcsprm *wcs, const double coord, |
---|
| 332 | std::string outputUnits) |
---|
[103] | 333 | { |
---|
[528] | 334 | /// @details |
---|
| 335 | /// Convert the wcs coordinate given by coord to a velocity in km/s. |
---|
| 336 | /// Does this by checking the ztype parameter in wcs to see if it is |
---|
| 337 | /// FREQ or otherwise, and converts accordingly. |
---|
| 338 | /// |
---|
| 339 | /// \param wcs The wcsprm struct containing the WCS information. |
---|
| 340 | /// \param coord The input WCS coordinate. |
---|
| 341 | /// \param outputUnits The string containing the units of the output velocity. |
---|
| 342 | /// \return The velocity of the input coord value. |
---|
[22] | 343 | |
---|
[120] | 344 | static int errflag = 0; |
---|
[103] | 345 | double scale, offset, power; |
---|
[160] | 346 | int specIndex = wcs->spec; |
---|
[478] | 347 | if(specIndex<0) specIndex = 2; |
---|
[492] | 348 | if(specIndex>=wcs->naxis) specIndex = wcs->naxis-1; |
---|
[205] | 349 | int status = wcsunits( wcs->cunit[specIndex], outputUnits.c_str(), |
---|
| 350 | &scale, &offset, &power); |
---|
[103] | 351 | |
---|
[205] | 352 | if(status > 0){ |
---|
[120] | 353 | if(errflag==0){ |
---|
| 354 | std::stringstream errmsg; |
---|
[205] | 355 | errmsg << "WCSUNITS Error Code = " << status << ":" |
---|
| 356 | << wcsunits_errmsg[status]; |
---|
| 357 | if(status == 10) errmsg << "\nTried to get conversion from " |
---|
| 358 | << wcs->cunit[specIndex] |
---|
| 359 | << " to " << outputUnits.c_str(); |
---|
[913] | 360 | // errmsg << "\nUsing coordinate value instead.\n"; |
---|
[1427] | 361 | DUCHAMPERROR("coordToVel", errmsg.str()); |
---|
[120] | 362 | errflag = 1; |
---|
[103] | 363 | } |
---|
| 364 | return coord; |
---|
[3] | 365 | } |
---|
[103] | 366 | else return pow( (coord*scale + offset), power); |
---|
| 367 | |
---|
[3] | 368 | } |
---|
| 369 | |
---|
[301] | 370 | double velToCoord(struct wcsprm *wcs, const float velocity, |
---|
| 371 | std::string outputUnits) |
---|
[3] | 372 | { |
---|
[528] | 373 | /// @details |
---|
| 374 | /// Convert the velocity given to the appropriate world coordinate for wcs. |
---|
| 375 | /// Does this by checking the ztype parameter in wcs to see if it is |
---|
| 376 | /// FREQ or otherwise, and converts accordingly. |
---|
| 377 | /// |
---|
| 378 | /// \param wcs The wcsprm struct containing the WCS information. |
---|
| 379 | /// \param velocity The input velocity. |
---|
| 380 | /// \param outputUnits The string containing the units of the input velocity. |
---|
| 381 | /// \return The WCS coordinate corresponding to the input velocity. |
---|
[22] | 382 | |
---|
[103] | 383 | double scale, offset, power; |
---|
[160] | 384 | int specIndex = wcs->spec; |
---|
[478] | 385 | if(specIndex<0) specIndex = 2; |
---|
[492] | 386 | if(specIndex>=wcs->naxis) specIndex = wcs->naxis-1; |
---|
[301] | 387 | int status = wcsunits( outputUnits.c_str(), wcs->cunit[specIndex], |
---|
[205] | 388 | &scale, &offset, &power); |
---|
[103] | 389 | |
---|
[205] | 390 | if(status > 0){ |
---|
[120] | 391 | std::stringstream errmsg; |
---|
[913] | 392 | errmsg << "WCSUNITS Error Code = " << status << ":" << wcsunits_errmsg[status]; |
---|
| 393 | // << "\nUsing coordinate value instead.\n"; |
---|
| 394 | DUCHAMPERROR("velToCoord",errmsg.str()); |
---|
[103] | 395 | return velocity; |
---|
| 396 | } |
---|
| 397 | else return (pow(velocity, 1./power) - offset) / scale; |
---|
[3] | 398 | |
---|
| 399 | } |
---|