diff options
-rw-r--r-- | tksao/frame/base.C | 7 | ||||
-rw-r--r-- | tksao/frame/fitsimage.C | 109 | ||||
-rw-r--r-- | tksao/frame/fitsimage.h | 6 | ||||
-rw-r--r-- | tksao/frame/fitsmap.C | 6 | ||||
-rw-r--r-- | tksao/vector/vector.C | 4 |
5 files changed, 72 insertions, 60 deletions
diff --git a/tksao/frame/base.C b/tksao/frame/base.C index 7667765..8f74b4b 100644 --- a/tksao/frame/base.C +++ b/tksao/frame/base.C @@ -642,7 +642,6 @@ Matrix Base::calcAlignWCS(FitsImage* fits1, FitsImage* fits2, Coord::CoordSystem sys1, Coord::CoordSystem sys2, Coord::SkyFrame sky) { - cerr << '*'; if ((!fits1 || !fits2) || (fits1 == fits2) || !(fits1->hasWCS(sys1)) || !(fits2->hasWCS(sys2))) return Matrix(); @@ -661,12 +660,10 @@ Matrix Base::calcAlignWCS(FitsImage* fits1, FitsImage* fits2, Matrix rr; if (cvt != AST__NULL) { astInvert(cvt); - Vector ll; Vector cc1 = fits1->center(); - fits1->astWCSTran(fits1->newast_, 1, cc1.v, cc1.v+1, 1, ll.v, ll.v+1); - Vector ur; + Vector ll = fits1->wcsTran(fits1->newast_, 1, cc1, 1); Vector cc2 = fits2->center(); - fits2->astWCSTran(fits2->newast_, 1, cc2.v, cc2.v+1, 1, ur.v, ur.v+1); + Vector ur = fits2->wcsTran(fits2->newast_, 1, cc2, 1); double fit[6]; double tol = 1; diff --git a/tksao/frame/fitsimage.C b/tksao/frame/fitsimage.C index 8805e7e..11fb0c2 100644 --- a/tksao/frame/fitsimage.C +++ b/tksao/frame/fitsimage.C @@ -1474,7 +1474,7 @@ void FitsImage::match(const char* xxname1, const char* yyname1, ixx1[ii] *= M_PI/180.; iyy1[ii] *= M_PI/180.; } - astWCSTran(ast_[ss], nxx1, ixx1, iyy1, 0, oxx1, oyy1); + wcsTran(ast_[ss], nxx1, ixx1, iyy1, 0, oxx1, oyy1); } } @@ -1489,7 +1489,7 @@ void FitsImage::match(const char* xxname1, const char* yyname1, ixx2[ii] *= M_PI/180.; iyy2[ii] *= M_PI/180.; } - astWCSTran(ast_[ss], nxx2, ixx2, iyy2, 0, oxx2, oyy2); + wcsTran(ast_[ss], nxx2, ixx2, iyy2, 0, oxx2, oyy2); } } @@ -1648,7 +1648,7 @@ void FitsImage::match(const char* xxname1, const char* yyname1, AstFrameSet* cvt = (AstFrameSet*)astConvert(wcs1, wcs2, "SKY"); if (cvt != AST__NULL) - astWCSTran(cvt, nxx1, ixx1, iyy1, 1, xx1, yy1); + wcsTran(cvt, nxx1, ixx1, iyy1, 1, xx1, yy1); } else { memcpy(xx1,ixx1,nxx1*sizeof(double)); @@ -2913,7 +2913,9 @@ Vector FitsImage::getWCScdelt(Coord::CoordSystem sys) else return Vector(); } -#else +#endif + +#ifdef NEWWCS double FitsImage::getWCSPixelSize(Coord::CoordSystem sys) { if (!hasWCS(sys)) @@ -2931,7 +2933,7 @@ double FitsImage::getWCSPixelSize(Coord::CoordSystem sys) yy[0] = cc[1]; yy[1] = cc[1]+1; yy[2] = cc[1]+1; - astWCSTran(newast_,3,xx,yy,1,wxx,wyy); + wcsTran(newast_,3,xx,yy,1,wxx,wyy); double pt0[2]; pt0[0] = wxx[0]; @@ -2967,7 +2969,7 @@ double FitsImage::getWCSPixelArea(Coord::CoordSystem sys) yy[0] = cc[1]; yy[1] = cc[1]+1; yy[2] = cc[1]+1; - astWCSTran(newast_,3,xx,yy,1,wxx,wyy); + wcsTran(newast_,3,xx,yy,1,wxx,wyy); double pt0[2]; pt0[0] = wxx[0]; @@ -3038,7 +3040,7 @@ Coord::Orientation FitsImage::getWCSOrientation(Coord::CoordSystem sys, yy[0] = pp[1]; yy[1] = pp[1]+1; yy[2] = pp[1]; - astWCSTran(newast_,3,xx,yy,1,wx,wy); + wcsTran(newast_,3,xx,yy,1,wx,wy); double aa[2], bb[2], cc[2]; aa[0]= wx[0]; @@ -3100,7 +3102,7 @@ double FitsImage::getWCSRotation(Coord::CoordSystem sys, Coord::SkyFrame sky) xx[1] = pp[0]; yy[0] = pp[1]; yy[1] = pp[1]+1; - astWCSTran(newast_,2,xx,yy,1,wx,wy); + wcsTran(newast_,2,xx,yy,1,wx,wy); double aa[2], bb[2]; aa[0]= wx[0]; @@ -3133,27 +3135,23 @@ const char* FitsImage::getWCSName(Coord::CoordSystem sys) #endif #ifndef NEWWCS -Vector FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys, +Vector FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys, Coord::SkyFrame sky) { astClearStatus; int ss = sys-Coord::WCS; if (ss>=0 && ast_ && ast_[ss]) { - double xx =0; - double yy =0; if (astWCSIsASkyFrame(ast_[ss])) { setAstWCSSkyFrame(ast_[ss],sky); - astWCSTran(ast_[ss], 1, in.v, in.v+1, 1, &xx, &yy); - if (astOK) - if (checkAstWCS(xx,yy)) - return Vector(radToDeg(xx),yy*180./M_PI); + Vector out = wcsTran(ast_[ss], in, 1); + if (astOK && checkWCS(out)) + return out.radToDeg(); } else { - astWCSTran(ast_[ss], 1, in.v, in.v+1, 1, &xx, &yy); - if (astOK) - if (checkAstWCS(xx,yy)) - return Vector(xx,yy); + Vector out = wcsTran(ast_[ss], in, 1); + if (astOK && checkWCS(out)) + return out; } } @@ -3176,8 +3174,8 @@ Vector FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys, double xx =0; double yy =0; - astWCSTran(newast_, 1, in.v, in.v+1, 1, &xx, &yy); - if (astOK && checkAstWCS(xx,yy)) { + wcsTran(newast_, 1, in.v, in.v+1, 1, &xx, &yy); + if (astOK && checkWCS(xx,yy)) { if (astWCSIsASkyFrame(newast_)) return Vector(radToDeg(xx),yy*180./M_PI); else @@ -3209,19 +3207,19 @@ Vector* FitsImage::pix2wcs(Vector* in, int num, Coord::CoordSystem sys, if (ss>=0 && ast_ && ast_[ss]) { if (astWCSIsASkyFrame(ast_[ss])) { setAstWCSSkyFrame(ast_[ss],sky); - astWCSTran(ast_[ss], num, xin, yin, 1, xout, yout); + wcsTran(ast_[ss], num, xin, yin, 1, xout, yout); if (astOK) { for (int ii=0; ii<num; ii++) - if (checkAstWCS(xout[ii],yout[ii])) + if (checkWCS(xout[ii],yout[ii])) out[ii] = Vector(radToDeg(xout[ii]),yout[ii]*180./M_PI); return out; } } else { - astWCSTran(ast_[ss], num, xin, yin, 1, xout, yout); + wcsTran(ast_[ss], num, xin, yin, 1, xout, yout); if (astOK) { for (int ii=0; ii<num; ii++) - if (checkAstWCS(xout[ii],yout[ii])) + if (checkWCS(xout[ii],yout[ii])) out[ii] = Vector(xout[ii],yout[ii]); return out; } @@ -3256,17 +3254,17 @@ Vector* FitsImage::pix2wcs(Vector* in, int num, Coord::CoordSystem sys, yin[ii] = (in[ii])[1]; } - astWCSTran(newast_, num, xin, yin, 1, xout, yout); + wcsTran(newast_, num, xin, yin, 1, xout, yout); if (astOK) { if (astWCSIsASkyFrame(newast_)) { for (int ii=0; ii<num; ii++) - if (checkAstWCS(xout[ii],yout[ii])) + if (checkWCS(xout[ii],yout[ii])) out[ii] = Vector(radToDeg(xout[ii]),yout[ii]*180./M_PI); return out; } else { for (int ii=0; ii<num; ii++) - if (checkAstWCS(xout[ii],yout[ii])) + if (checkWCS(xout[ii],yout[ii])) out[ii] = Vector(xout[ii],yout[ii]); return out; } @@ -3291,8 +3289,8 @@ char* FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys, ostringstream str; if (astWCSIsASkyFrame(ast_[ss])) { setAstWCSSkyFrame(ast_[ss],sky); - astWCSTran(ast_[ss], 1, in.v, in.v+1, 1, &xx, &yy); - if (!astOK || !checkAstWCS(xx,yy)) { + wcsTran(ast_[ss], 1, in.v, in.v+1, 1, &xx, &yy); + if (!astOK || !checkWCS(xx,yy)) { maperr =1; lbuf[0] = '\0'; return lbuf; @@ -3333,8 +3331,8 @@ char* FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys, } } else { - astWCSTran(ast_[ss], 1, in.v, in.v+1, 1, &xx, &yy); - if (!astOK || !checkAstWCS(xx,yy)) { + wcsTran(ast_[ss], 1, in.v, in.v+1, 1, &xx, &yy); + if (!astOK || !checkWCS(xx,yy)) { maperr =1; return lbuf; } @@ -3368,8 +3366,8 @@ char* FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys, double yy =0; ostringstream str; - astWCSTran(newast_, 1, in.v, in.v+1, 1, &xx, &yy); - if (astOK && checkAstWCS(xx,yy)) { + wcsTran(newast_, 1, in.v, in.v+1, 1, &xx, &yy); + if (astOK && checkWCS(xx,yy)) { if (astWCSIsASkyFrame(newast_)) { switch (format) { case Coord::DEGREES: @@ -3431,15 +3429,15 @@ Vector FitsImage::wcs2pix(Vector in, Coord::CoordSystem sys, if (astWCSIsASkyFrame(ast_[ss])) { setAstWCSSkyFrame(ast_[ss],sky); Vector rr = in*M_PI/180.; - astWCSTran(ast_[ss], 1, rr.v, &(rr[1]), 0, &xx, &yy); + wcsTran(ast_[ss], 1, rr.v, &(rr[1]), 0, &xx, &yy); if (astOK) - if (checkAstWCS(xx,yy)) + if (checkWCS(xx,yy)) return Vector(xx,yy); } else { - astWCSTran(ast_[ss], 1, in.v, in.v+1, 0, &xx, &yy); + wcsTran(ast_[ss], 1, in.v, in.v+1, 0, &xx, &yy); if (astOK) - if (checkAstWCS(xx,yy)) + if (checkWCS(xx,yy)) return Vector(xx,yy); } } @@ -3467,8 +3465,8 @@ Vector FitsImage::wcs2pix(Vector in, Coord::CoordSystem sys, if (astWCSIsASkyFrame(newast_)) in *= M_PI/180.; - astWCSTran(newast_, 1, in.v, in.v+1, 0, &xx, &yy); - if (astOK && checkAstWCS(xx,yy)) + wcsTran(newast_, 1, in.v, in.v+1, 0, &xx, &yy); + if (astOK && checkWCS(xx,yy)) return Vector(xx,yy); maperr =1; @@ -3500,19 +3498,19 @@ Vector* FitsImage::wcs2pix(Vector* in, int num, Coord::CoordSystem sys, xin[kk] *= M_PI/180.; yin[kk] *= M_PI/180.; } - astWCSTran(ast_[ss], num, xin, yin, 0, xout, yout); + wcsTran(ast_[ss], num, xin, yin, 0, xout, yout); if (astOK) { for (int kk=0; kk<num; kk++) - if (checkAstWCS(xout[kk],yout[kk])) + if (checkWCS(xout[kk],yout[kk])) out[kk] = Vector(xout[kk],yout[kk]); return out; } } else { - astWCSTran(ast_[ss], num, xin, yin, 0, xout, yout); + wcsTran(ast_[ss], num, xin, yin, 0, xout, yout); if (astOK) { for (int kk=0; kk<num; kk++) - if (checkAstWCS(xout[kk],yout[kk])) + if (checkWCS(xout[kk],yout[kk])) out[kk] = Vector(xout[kk],yout[kk]); return out; } @@ -3554,10 +3552,10 @@ Vector* FitsImage::wcs2pix(Vector* in, int num, Coord::CoordSystem sys, } } - astWCSTran(newast_, num, xin, yin, 0, xout, yout); + wcsTran(newast_, num, xin, yin, 0, xout, yout); if (astOK) { for (int kk=0; kk<num; kk++) - if (checkAstWCS(xout[kk],yout[kk])) + if (checkWCS(xout[kk],yout[kk])) out[kk] = Vector(xout[kk],yout[kk]); return out; } @@ -4050,12 +4048,18 @@ void FitsImage::astinit0(int ss, FitsHead* hd, FitsHead* prim) setAstWCSSkyFrame(ast_[ss],Coord::FK5); } -int FitsImage::checkAstWCS(double xx, double yy) +int FitsImage::checkWCS(double xx, double yy) { // check for reasonable values return (fabs(xx) < FLT_MAX && fabs(yy) < FLT_MAX) ? 1 : 0; } +int FitsImage::checkWCS(Vector& vv) +{ + // check for reasonable values + return (fabs(vv[0]) < FLT_MAX && fabs(vv[1]) < FLT_MAX) ? 1 : 0; +} + void FitsImage::setAstWCSFormat(AstFrameSet* aa, int id, const char* format) { // is it already set? @@ -4197,15 +4201,22 @@ int FitsImage::astWCSIsASkyFrame(void* ast) #endif #ifndef NEWWCS -void FitsImage::astWCSTran(AstFrameSet* ast, int npoint, +void FitsImage::wcsTran(AstFrameSet* ast, int npoint, const double* xin, const double* yin, int forward, double* xout, double* yout) { astTran2(ast, npoint, xin, yin, forward, xout, yout); } + +Vector FitsImage::wcsTran(AstFrameSet* ast, Vector& in, int forward) +{ + double xout, yout; + astTran2(ast, 1, in.v, in.v+1, forward, &xout, &yout); + return Vector(xout, yout); +} #else -void FitsImage::astWCSTran(AstFrameSet* ast, int npoint, +void FitsImage::wcsTran(AstFrameSet* ast, int npoint, const double* xin, const double* yin, int forward, double* xout, double* yout) diff --git a/tksao/frame/fitsimage.h b/tksao/frame/fitsimage.h index 3759285..2765681 100644 --- a/tksao/frame/fitsimage.h +++ b/tksao/frame/fitsimage.h @@ -148,7 +148,8 @@ class FitsImage { void astinit(FitsHead*, FitsHead*); #endif void astinit0(int, FitsHead*, FitsHead*); - int checkAstWCS(double, double); + int checkWCS(double, double); + int checkWCS(Vector&); AstFrameSet* fits2ast(FitsHead*); AstFrameSet* buildast(int, FitsHead*, FitsHead*); AstFrameSet* buildast0(int, FitsHead*, FitsHead*); @@ -401,8 +402,9 @@ class FitsImage { #endif int astWCSIsASkyFrame(void*); - void astWCSTran(AstFrameSet*, int npoint, const double*, const double*, + void wcsTran(AstFrameSet*, int npoint, const double*, const double*, int, double*, double*); + Vector wcsTran(AstFrameSet* ast, Vector& in, int forward); #ifdef NEWWCS void setAstWCSSystem(AstFrameSet*, Coord::CoordSystem); diff --git a/tksao/frame/fitsmap.C b/tksao/frame/fitsmap.C index 6f666f3..6f0e975 100644 --- a/tksao/frame/fitsmap.C +++ b/tksao/frame/fitsmap.C @@ -185,7 +185,7 @@ double FitsImage::mapLenFromRef(double dd, Coord::CoordSystem sys, double yy[2], wyy[2]; yy[0] = cc[1]; yy[1] = cc[1]+dd; - astWCSTran(newast_,2,xx,yy,1,wxx,wyy); + wcsTran(newast_,2,xx,yy,1,wxx,wyy); double pt0[2]; pt0[0] = wxx[0]; @@ -311,7 +311,7 @@ double FitsImage::mapLenToRef(double dd, Coord::CoordSystem sys, Vector cc = center(); Vector wcc; - astWCSTran(newast_,1,cc.v,cc.v+1,1,wcc.v,wcc.v+1); + wcsTran(newast_,1,cc.v,cc.v+1,1,wcc.v,wcc.v+1); double wxx[2], xx[2]; wxx[0] = wcc[0]; @@ -319,7 +319,7 @@ double FitsImage::mapLenToRef(double dd, Coord::CoordSystem sys, double wyy[2], yy[2]; wyy[0] = wcc[1]; wyy[1] = wcc[1]+rdd; - astWCSTran(newast_,2,wxx,wyy,0,xx,yy); + wcsTran(newast_,2,wxx,wyy,0,xx,yy); double pt0[2]; pt0[0] = xx[0]; diff --git a/tksao/vector/vector.C b/tksao/vector/vector.C index 65ca1c3..429fcc0 100644 --- a/tksao/vector/vector.C +++ b/tksao/vector/vector.C @@ -7,6 +7,7 @@ #include "vector.h" #include "vector3d.h" #include "fuzzy.h" +#include "util.h" // Vector::manip @@ -46,7 +47,8 @@ Vector& Vector::clip(const BBox& bb) Vector Vector::radToDeg() { - return Vector(v[0]*180./M_PI,v[1]*180./M_PI); + // we want the first coord to be 0-360 + return Vector(::radToDeg(v[0]),v[1]*180./M_PI); } Vector Vector::degToRad() |