From da35ec15dc1ff6bb6331f83ef186c1c96c1aa484 Mon Sep 17 00:00:00 2001 From: William Joye Date: Tue, 7 Aug 2018 17:24:44 -0400 Subject: simplify wcs code --- tksao/frame/base.C | 7 +- tksao/frame/fitsimage.C | 258 ++++++++++++++++++++++++++++++++++++------------ tksao/frame/fitsimage.h | 15 +-- tksao/frame/grid25d.C | 3 +- tksao/frame/grid2d.C | 3 +- tksao/frame/grid3d.C | 3 +- 6 files changed, 213 insertions(+), 76 deletions(-) diff --git a/tksao/frame/base.C b/tksao/frame/base.C index 43e2c8d..afacd6b 100644 --- a/tksao/frame/base.C +++ b/tksao/frame/base.C @@ -387,11 +387,14 @@ Matrix Base::calcAlignWCS(FitsImage* fits1, FitsImage* fits2, astClearStatus; // just to make sure astBegin; // start memory management - fits1->setWCSSkyFrame(sys1, sky); AstFrameSet* wcs1 = (AstFrameSet*)astCopy(fits1->ast_); + fits1->wcsSystem(wcs1,sys1); + fits1->wcsSkyFrame(wcs1,sky); astInvert(wcs1); - fits2->setWCSSkyFrame(sys2, sky); + AstFrameSet* wcs2 = (AstFrameSet*)astCopy(fits2->ast_); + fits2->wcsSystem(wcs2,sys1); + fits2->wcsSkyFrame(wcs2,sky); astInvert(wcs2); AstFrameSet* cvt = (AstFrameSet*)astConvert(wcs2, wcs1, ""); diff --git a/tksao/frame/fitsimage.C b/tksao/frame/fitsimage.C index 328ed55..e3d6a4a 100644 --- a/tksao/frame/fitsimage.C +++ b/tksao/frame/fitsimage.C @@ -1328,11 +1328,13 @@ void FitsImage::match(const char* xxname1, const char* yyname1, Vector* ptr1 =NULL; if (sky1 != sky2) { - setWCSSkyFrame(sys1, sky1); AstFrameSet* wcs1 = (AstFrameSet*)astCopy(ast_); + wcsSystem(wcs1,sys1); + wcsSkyFrame(wcs1,sky1); - setWCSSkyFrame(sys2, sky2); AstFrameSet* wcs2 = (AstFrameSet*)astCopy(ast_); + wcsSystem(wcs2,sys2); + wcsSkyFrame(wcs2,sky2); AstFrameSet* cvt = (AstFrameSet*)astConvert(wcs1, wcs2, "SKY"); if (cvt != AST__NULL) { @@ -1345,11 +1347,14 @@ void FitsImage::match(const char* xxname1, const char* yyname1, // now compare if (ptr1 && ptr2) { - setWCSSkyFrame(sys2, sky2); + AstFrameSet* wcs = (AstFrameSet*)astCopy(ast_); + wcsSystem(wcs,sys2); + wcsSkyFrame(wcs,sky2); + Tcl_Obj* objrr = Tcl_NewListObj(0,NULL); for(int jj=0; jj0) || (!hasWCSCel(sys) && ang<0)) @@ -2616,9 +2621,9 @@ double FitsImage::getWCSRotation(Coord::CoordSystem sys, Coord::SkyFrame sky) in[0] = center(); in[1] = center()+Vector(0,1); in[2] = center()+Vector(1,0); - wcsTran(3, in, 1, out); - double ang = wcsAxAngle(out[0], out[1]); - double ang3 = wcsAngle(out[1],out[0],out[2]); + wcsTran(ast_, 3, in, 1, out); + double ang = wcsAxAngle(ast_,out[0], out[1]); + double ang3 = wcsAngle(ast_,out[1],out[0],out[2]); if ((!(isnan(ang )||isinf(ang) ||(ang ==-DBL_MAX)||(ang ==DBL_MAX))) && (!(isnan(ang3)||isinf(ang3)||(ang3==-DBL_MAX)||(ang3==DBL_MAX)))) { @@ -2632,14 +2637,14 @@ double FitsImage::getWCSRotation(Coord::CoordSystem sys, Coord::SkyFrame sky) } else { // special case for HPX Vector cc = center(); - Vector wcc = wcsTran(cc, 1); + Vector wcc = wcsTran(ast_, cc, 1); Vector wnorth = wcc + Vector(0,.001); - Vector north = wcsTran(wnorth,0); + Vector north = wcsTran(ast_, wnorth,0); int current = astGetI(ast_,"Current"); int base = astGetI(ast_,"Base"); astSetI(ast_,"Current",base); - double ang = wcsAxAngle(cc,north); + double ang = wcsAxAngle(ast_,cc,north); astSetI(ast_,"Current",current); if (!(isnan(ang)||isinf(ang)||(ang == -DBL_MAX)||(ang == DBL_MAX))) return ang; @@ -2666,7 +2671,7 @@ Vector FitsImage::pix2wcs(const Vector& in, Coord::CoordSystem sys, setWCSSkyFrame(sys, sky); - Vector out = wcsTran(in, 1); + Vector out = wcsTran(ast_, in, 1); if (astOK && checkWCS(out)) return hasWCSCel(sys) ? zero360(radToDeg(out)) : out; else @@ -2686,7 +2691,7 @@ char* FitsImage::pix2wcs(const Vector& in, Coord::CoordSystem sys, setWCSSkyFrame(sys, sky); ostringstream str; - Vector out = wcsTran(in, 1); + Vector out = wcsTran(ast_, in, 1); if (astOK && checkWCS(out)) { if (hasWCSCel(sys)) { switch (format) { @@ -2744,7 +2749,7 @@ Vector3d FitsImage::pix2wcs(const Vector3d& in, Coord::CoordSystem sys, setWCSSkyFrame(sys, sky); - Vector3d out = wcsTran(in, 1); + Vector3d out = wcsTran(ast_, in, 1); if (astOK && checkWCS(out)) return hasWCSCel(sys) ? zero360(radToDeg(out)) : out; else @@ -2764,7 +2769,7 @@ char* FitsImage::pix2wcs(const Vector3d& in, Coord::CoordSystem sys, setWCSSkyFrame(sys, sky); ostringstream str; - Vector3d out = wcsTran(in, 1); + Vector3d out = wcsTran(ast_, in, 1); if (astOK && checkWCS(out)) { if (hasWCSCel(sys)) { switch (format) { @@ -2820,7 +2825,7 @@ Vector FitsImage::wcs2pix(const Vector& vv, Coord::CoordSystem sys, setWCSSkyFrame(sys, sky); Vector in = hasWCSCel(sys) ? degToRad(vv) : vv; - Vector out = wcsTran(in, 0); + Vector out = wcsTran(ast_, in, 0); if (astOK && checkWCS(out)) return out; } @@ -2838,7 +2843,7 @@ Vector3d FitsImage::wcs2pix(const Vector3d& vv, Coord::CoordSystem sys, setWCSSkyFrame(sys, sky); Vector3d in = hasWCSCel(sys) ? degToRad(vv) : vv; - Vector3d out = wcsTran(in, 0); + Vector3d out = wcsTran(ast_, in, 0); if (astOK && checkWCS(out)) return out; } @@ -2856,8 +2861,8 @@ double FitsImage::getWCSDist(const Vector& vv1, const Vector& vv2, setWCSSkyFrame(sys, Coord::FK5); return hasWCSCel(sys) ? - radToDeg(wcsDistance(degToRad(vv1), degToRad(vv2))) : - wcsDistance(vv1, vv2); + radToDeg(wcsDistance(ast_,degToRad(vv1),degToRad(vv2))) : + wcsDistance(ast_,vv1,vv2); } int FitsImage::hasWCS(Coord::CoordSystem sys) @@ -3096,8 +3101,77 @@ int FitsImage::checkWCS(Vector3d& vv) fabs(vv[2]) < FLT_MAX ) ? 1 : 0; } +#if 0 +void FitsImage::setASTFormat(Coord::CoordSystem sys, Coord::SkyFrame sky, + Coord::SkyFormat format) +{ + int nn = astGetI(ast_, "Nframe"); + // do we have a AST wcs? + if (hasWCSAST) { + for (int ii=0; iiparent_->precDeg_; + setWCSFormat(1,str.str().c_str()); + setWCSFormat(2,str.str().c_str()); + break; + + case Coord::SEXAGESIMAL: + { + hms << "hms." << context_->parent_->precHMS_; + dms << "+dms." << context_->parent_->precDMS_; + switch (sky) { + case Coord::FK4: + case Coord::FK5: + case Coord::ICRS: + setWCSFormat(1,hms.str().c_str()); + setWCSFormat(2,dms.str().c_str()); + break; + case Coord::GALACTIC: + case Coord::ECLIPTIC: + setWCSFormat(1,dms.str().c_str()); + setWCSFormat(2,dms.str().c_str()); + break; + } + } + break; + } + } + } + } + else { + str << "1." << context_->parent_->precLinear_ << 'G'; + setWCSFormat(1,str.str().c_str()); + setWCSFormat(2,str.str().c_str()); + } + } +} +#endif + void FitsImage::setWCSFormat(int id, const char* format) { + // cerr << 'f'; // is it already set? // ast is very slow when changing params { @@ -3115,6 +3189,7 @@ void FitsImage::setWCSFormat(int id, const char* format) void FitsImage::setWCSSkyFrame(Coord::CoordSystem sys, Coord::SkyFrame sky) { + // cerr << 's'; int nn = astGetI(ast_,"nframe"); char cc = ' '; int ww = sys-Coord::WCS; @@ -3179,41 +3254,70 @@ void FitsImage::setWCSSkyFrame(Coord::CoordSystem sys, Coord::SkyFrame sky) astSetD(ast_, "EQUINOX", astGetD(ast_, "EPOCH")); break; } +} -#if 0 - for (int ii=0; iislice(2) : 0; - astTranN(ast_, 1, 3, 1, pin, forward, 3, 1, pout); + astTranN(ast, 1, 3, 1, pin, forward, 3, 1, pout); return Vector(pout[0],pout[1]); } break; @@ -3234,7 +3338,7 @@ Vector FitsImage::wcsTran(const Vector& in, int forward) pin[1] = in[1]; pin[2] = forward ? context_->slice(2) : 0; pin[3] = forward ? context_->slice(3) : 0; - astTranN(ast_, 1, 4, 1, pin, forward, 4, 1, pout); + astTranN(ast, 1, 4, 1, pin, forward, 4, 1, pout); return Vector(pout[0],pout[1]); } break; @@ -3355,9 +3459,9 @@ void FitsImage::wcsTran(AstFrameSet* ast, int npoint, } } -Vector3d FitsImage::wcsTran(const Vector3d& in, int forward) +Vector3d FitsImage::wcsTran(AstFrameSet* ast, const Vector3d& in, int forward) { - int naxes = astGetI(ast_,"Naxes"); + int naxes = astGetI(ast,"Naxes"); switch (naxes) { case 1: case 2: @@ -3365,7 +3469,7 @@ Vector3d FitsImage::wcsTran(const Vector3d& in, int forward) double pout[2]; pin[0] = in[0]; pin[1] = in[1]; - astTranN(ast_, 1, 2, 1, pin, forward, 2, 1, pout); + astTranN(ast, 1, 2, 1, pin, forward, 2, 1, pout); return Vector3d(pout[0],pout[1],forward ? 1 : 0); break; case 3: @@ -3375,7 +3479,7 @@ Vector3d FitsImage::wcsTran(const Vector3d& in, int forward) pin[0] = in[0]; pin[1] = in[1]; pin[2] = in[2]; - astTranN(ast_, 1, 3, 1, pin, forward, 3, 1, pout); + astTranN(ast, 1, 3, 1, pin, forward, 3, 1, pout); return Vector3d(pout[0],pout[1],pout[2]); } break; @@ -3387,7 +3491,7 @@ Vector3d FitsImage::wcsTran(const Vector3d& in, int forward) pin[1] = in[1]; pin[2] = in[2]; pin[3] = forward ? context_->slice(3) : 0; - astTranN(ast_, 1, 4, 1, pin, forward, 4, 1, pout); + astTranN(ast, 1, 4, 1, pin, forward, 4, 1, pout); return Vector3d(pout[0],pout[1],pout[2]); } break; @@ -3395,15 +3499,16 @@ Vector3d FitsImage::wcsTran(const Vector3d& in, int forward) return Vector3d(); } -double FitsImage::wcsDistance(const Vector& vv1, const Vector& vv2) +double FitsImage::wcsDistance(AstFrameSet* ast, + const Vector& vv1, const Vector& vv2) { - int naxes = astGetI(ast_,"Naxes"); + int naxes = astGetI(ast,"Naxes"); switch (naxes) { case 1: // error break; case 2: - return astDistance(ast_, vv1.v, vv2.v); + return astDistance(ast, vv1.v, vv2.v); case 3: { double ptr1[3]; @@ -3415,7 +3520,7 @@ double FitsImage::wcsDistance(const Vector& vv1, const Vector& vv2) ptr2[1] = vv2[1]; ptr2[2] = 0; - return astDistance(ast_, ptr1, ptr2); + return astDistance(ast, ptr1, ptr2); } case 4: { @@ -3430,23 +3535,24 @@ double FitsImage::wcsDistance(const Vector& vv1, const Vector& vv2) ptr2[2] = 0; ptr2[3] = 0; - return astDistance(ast_, ptr1, ptr2); + return astDistance(ast, ptr1, ptr2); } } return 0; } -double FitsImage::wcsAngle(const Vector& vv1, const Vector& vv2, +double FitsImage::wcsAngle(AstFrameSet* ast, + const Vector& vv1, const Vector& vv2, const Vector& vv3) { - int naxes = astGetI(ast_,"Naxes"); + int naxes = astGetI(ast,"Naxes"); switch (naxes) { case 1: // error break; case 2: - return astAngle(ast_,vv1.v,vv2.v,vv3.v); + return astAngle(ast,vv1.v,vv2.v,vv3.v); case 3: { double ptr1[3]; @@ -3462,7 +3568,7 @@ double FitsImage::wcsAngle(const Vector& vv1, const Vector& vv2, ptr3[1] = vv3[1]; ptr3[2] = 0; - return astAngle(ast_, ptr1, ptr2, ptr3); + return astAngle(ast, ptr1, ptr2, ptr3); } case 4: { @@ -3482,22 +3588,23 @@ double FitsImage::wcsAngle(const Vector& vv1, const Vector& vv2, ptr3[2] = 0; ptr3[3] = 0; - return astAngle(ast_, ptr1, ptr2, ptr3); + return astAngle(ast, ptr1, ptr2, ptr3); } } return 0; } -double FitsImage::wcsAxAngle(const Vector& vv1, const Vector& vv2) +double FitsImage::wcsAxAngle(AstFrameSet* ast, + const Vector& vv1, const Vector& vv2) { - int naxes = astGetI(ast_,"Naxes"); + int naxes = astGetI(ast,"Naxes"); switch (naxes) { case 1: // error break; case 2: - return astAxAngle(ast_, vv1.v, vv2.v, 2); + return astAxAngle(ast, vv1.v, vv2.v, 2); case 3: { double ptr1[3]; @@ -3509,7 +3616,7 @@ double FitsImage::wcsAxAngle(const Vector& vv1, const Vector& vv2) ptr2[1] = vv2[1]; ptr2[2] = 0; - return astAxAngle(ast_, ptr1, ptr2, 2); + return astAxAngle(ast, ptr1, ptr2, 2); } case 4: { @@ -3524,7 +3631,7 @@ double FitsImage::wcsAxAngle(const Vector& vv1, const Vector& vv2) ptr2[2] = 0; ptr2[3] = 0; - return astAxAngle(ast_, ptr1, ptr2, 2); + return astAxAngle(ast, ptr1, ptr2, 2); } } @@ -3706,3 +3813,26 @@ AstFrameSet* FitsImage::fits2ast(FitsHead* hd) return frameSet; } + +#if 0 +void FitsImage::dump() +{ + AstFrameSet* fs = + (AstFrameSet*)astFindFrame(ast_, astSkyFrame(" MaxAxes=4")," "); + if (fs) { + for (int ii=0; ii<2; ii++) { + ostringstream str3; + str3 << "Format(" << ii+1 << ")" << ends; + + ostringstream str4; + str4 << "Symbol(" << ii+1 << ")" << ends; + + cerr << ii << ' ' + << astGetC(ast_, str3.str().c_str()) << ' ' + << astGetC(ast_, str4.str().c_str()) << ' ' + << endl; + } + cerr << astGetC(fs, "System") << ' ' << endl; + } +} +#endif diff --git a/tksao/frame/fitsimage.h b/tksao/frame/fitsimage.h index 4f953b6..e44ee31 100644 --- a/tksao/frame/fitsimage.h +++ b/tksao/frame/fitsimage.h @@ -389,14 +389,15 @@ class FitsImage { double getWCSDist(const Vector&, const Vector&, Coord::CoordSystem); const char* getWCSName(Coord::CoordSystem); - Vector wcsTran(const Vector&, int); - Vector3d wcsTran(const Vector3d&, int); - void wcsTran(int num, Vector* in, int forward, Vector* out) - {wcsTran(ast_,num,in,forward,out);} + Vector wcsTran(AstFrameSet*, const Vector&, int); void wcsTran(AstFrameSet*, int, Vector*, int, Vector*); - double wcsDistance(const Vector&, const Vector&); - double wcsAngle(const Vector&, const Vector&, const Vector&); - double wcsAxAngle(const Vector&, const Vector&); + Vector3d wcsTran(AstFrameSet*, const Vector3d&, int); + double wcsDistance(AstFrameSet*, const Vector&, const Vector&); + double wcsAngle(AstFrameSet*, const Vector&, const Vector&, const Vector&); + double wcsAxAngle(AstFrameSet*, const Vector&, const Vector&); + + int wcsSystem(AstFrameSet* ast, Coord::CoordSystem sys); + int wcsSkyFrame(AstFrameSet* ast, Coord::SkyFrame sky); void setWCSSkyFrame(Coord::CoordSystem, Coord::SkyFrame); void setWCSFormat(int, const char*); diff --git a/tksao/frame/grid25d.C b/tksao/frame/grid25d.C index ae5a428..17956e7 100644 --- a/tksao/frame/grid25d.C +++ b/tksao/frame/grid25d.C @@ -63,8 +63,9 @@ int Grid25d::doit(RenderMode rm) return 1; } - fits->setWCSSkyFrame(system_, sky_); AstFrameSet* ast = (AstFrameSet*)astCopy(fits->ast_); + fits->wcsSystem(ast,system_); + fits->wcsSkyFrame(ast,sky_); int naxes = astGetI(ast,"Naxes"); switch (naxes) { diff --git a/tksao/frame/grid2d.C b/tksao/frame/grid2d.C index 567a3f8..02a80c7 100644 --- a/tksao/frame/grid2d.C +++ b/tksao/frame/grid2d.C @@ -64,8 +64,9 @@ int Grid2d::doit(RenderMode rm) return 1; } - fits->setWCSSkyFrame(system_, sky_); AstFrameSet* ast = (AstFrameSet*)astCopy(fits->ast_); + fits->wcsSystem(ast,system_); + fits->wcsSkyFrame(ast,sky_); int naxes = astGetI(ast,"Naxes"); switch (naxes) { diff --git a/tksao/frame/grid3d.C b/tksao/frame/grid3d.C index e36b36d..8dad14b 100644 --- a/tksao/frame/grid3d.C +++ b/tksao/frame/grid3d.C @@ -64,8 +64,9 @@ int Grid3d::doit(RenderMode rm) return 1; } - fits->setWCSSkyFrame(system_, sky_); AstFrameSet* ast = (AstFrameSet*)astCopy(fits->ast_); + fits->wcsSystem(ast,system_); + fits->wcsSkyFrame(ast,sky_); int naxes = astGetI(ast,"Naxes"); switch (naxes) { -- cgit v0.12