diff options
| -rw-r--r-- | tksao/frame/base.C | 67 | ||||
| -rw-r--r-- | tksao/frame/fitsimage.C | 229 | ||||
| -rw-r--r-- | tksao/frame/fitsimage.h | 8 |
3 files changed, 240 insertions, 64 deletions
diff --git a/tksao/frame/base.C b/tksao/frame/base.C index bd54ed6..1c87af8 100644 --- a/tksao/frame/base.C +++ b/tksao/frame/base.C @@ -642,48 +642,45 @@ Matrix Base::calcAlignWCS(FitsImage* fits1, FitsImage* fits2, Coord::CoordSystem sys1, Coord::CoordSystem sys2, Coord::SkyFrame sky) { - if ((!fits1 || !fits2) || - (fits1 == fits2) || - !(fits1->hasWCS(sys1)) || - !(fits2->hasWCS(sys2))) + cerr << '*'; + if ((!fits1 || !fits2) || (fits1 == fits2) || + !(fits1->hasWCS(sys1)) || !(fits2->hasWCS(sys2))) return Matrix(); astClearStatus; // just to make sure astBegin; // start memory management - int ss1 = sys1-Coord::WCS; - int ss2 = sys2-Coord::WCS; + fits1->setAstWCSSystem(fits1->newast_, sys1); + fits2->setAstWCSSystem(fits2->newast_, sys2); + + AstFrameSet* wcs1 = (AstFrameSet*)astCopy(fits1->newast_); + AstFrameSet* wcs2 = (AstFrameSet*)astCopy(fits2->newast_); + astInvert(wcs1); + astInvert(wcs2); + AstFrameSet* cvt = (AstFrameSet*)astConvert(wcs1, wcs2, ""); Matrix rr; - if ((ss1>=0 && fits1->ast_ && fits1->ast_[ss1]) && - (ss2>=0 && fits2->ast_ && fits2->ast_[ss2])) { - AstFrameSet* wcs1 = (AstFrameSet*)astCopy(fits1->ast_[ss1]); - AstFrameSet* wcs2 = (AstFrameSet*)astCopy(fits2->ast_[ss2]); - astInvert(wcs1); - astInvert(wcs2); - AstFrameSet* cvt = (AstFrameSet*)astConvert(wcs1, wcs2, ""); - if (cvt != AST__NULL) { - astInvert(cvt); - Vector ll; - Vector cc1 = fits1->center(); - fits1->astWCSTran(fits1->ast_[ss1], 1, cc1.v, cc1.v+1, 1, ll.v, ll.v+1); - Vector ur; - Vector cc2 = fits2->center(); - fits2->astWCSTran(fits2->ast_[ss2], 1, cc2.v, cc2.v+1, 1, ur.v, ur.v+1); - - double fit[6]; - double tol = 1; - if (astLinearApprox(cvt, ll.v, ur.v, tol, fit) != AST__BAD) { - // fix the fit from AST - if (fit[2] == 0) { - fit[2] =1; - fit[0] =0; - } - if (fit[5] ==0) { - fit[5] =1; - fit[1] =0; - } - rr = Matrix(fit[2],fit[4],fit[3],fit[5],fit[0],fit[1]); + 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 cc2 = fits2->center(); + fits2->astWCSTran(fits2->newast_, 1, cc2.v, cc2.v+1, 1, ur.v, ur.v+1); + + double fit[6]; + double tol = 1; + if (astLinearApprox(cvt, ll.v, ur.v, tol, fit) != AST__BAD) { + // fix the fit from AST + if (fit[2] == 0) { + fit[2] =1; + fit[0] =0; + } + if (fit[5] ==0) { + fit[5] =1; + fit[1] =0; } + rr = Matrix(fit[2],fit[4],fit[3],fit[5],fit[0],fit[1]); } } diff --git a/tksao/frame/fitsimage.C b/tksao/frame/fitsimage.C index c349dc2..6282a85 100644 --- a/tksao/frame/fitsimage.C +++ b/tksao/frame/fitsimage.C @@ -109,6 +109,7 @@ FitsImage::FitsImage(Context* cx, Tcl_Interp* pp) wcs_ =NULL; wcsx_ =NULL; ast_ =NULL; + newast_ =NULL; wcsHeader_ =NULL; altHeader_ =NULL; @@ -184,6 +185,10 @@ FitsImage::~FitsImage() astAnnul(ast_[ii]); delete [] ast_; } +#ifdef NEWWCS + if (manageWCS_ && newast_) + astAnnul(newast_); +#endif if (wcsHeader_) delete wcsHeader_; @@ -1080,6 +1085,11 @@ void FitsImage::initWCS() ast_ = new AstFrameSet*[MULTWCSA]; for (int ii=0; ii<MULTWCSA; ii++) ast_[ii] = NULL; +#ifdef NEWWCS + if (manageWCS_ && newast_) + astAnnul(newast_); + newast_ = NULL; +#endif // shareWCS? manageWCS_ =1; @@ -1098,6 +1108,9 @@ void FitsImage::initWCS() wcsx_[ii] = ptr->wcsx_[ii]; for (int ii=0; ii<MULTWCSA; ii++) ast_[ii] = ptr->ast_[ii]; +#ifdef NEWWCS + newast_ = ptr->newast_; +#endif #ifndef NEWWCS initWCSPhysical(); @@ -1165,10 +1178,17 @@ void FitsImage::initWCS() astinit(ii, hd, prim); +#ifndef NEWWCS if (DebugAST) astShow(ast_[ii]); +#endif } } +#ifdef NEWWCS + astinit(hd, prim); + if (DebugAST && newast_) + astShow(newast_); +#endif #ifndef NEWWCS // WCSDEP @@ -1529,7 +1549,7 @@ void FitsImage::match(const char* xxname1, const char* yyname1, const char* rrname) { // only good for skyframe - + astClearStatus; // just to make sure astBegin; // start memory management @@ -1581,16 +1601,13 @@ void FitsImage::match(const char* xxname1, const char* yyname1, for (int ii=0 ; ii<nyy2 ; ii++) Tcl_GetDoubleFromObj(interp_, objyy2[ii], yy2+ii); - if (!hasWCS(sys1)) - return; - if (!hasWCS(sys2)) + if (!hasWCS(sys1) || !hasWCS(sys2)) return; int ss1 = sys1-Coord::WCS; int ss2 = sys2-Coord::WCS; // are both skyframe? - if (!((astWCSIsASkyFrame(ast_[ss1]) && - (astWCSIsASkyFrame(ast_[ss2]))))) + if (!((astWCSIsASkyFrame(ast_[ss1]) && (astWCSIsASkyFrame(ast_[ss2]))))) return; setAstWCSSkyFrame(ast_[ss1],sky1); @@ -3351,13 +3368,14 @@ double FitsImage::getWCSDist(Vector a, Vector b, Coord::CoordSystem sys) return rr; } +#ifndef NEWWCS + int FitsImage::hasWCS(Coord::CoordSystem sys) { int ss = sys-Coord::WCS; return (sys>=Coord::WCS && ast_ && ast_[ss]) ? 1 : 0; } -#ifndef NEWWCS int FitsImage::hasWCSEqu(Coord::CoordSystem sys) { astClearStatus; @@ -3389,39 +3407,146 @@ int FitsImage::hasWCSEqu(Coord::CoordSystem sys) return 0; } + +int FitsImage::hasWCSCel(Coord::CoordSystem sys) +{ + astClearStatus; + + int ss = sys-Coord::WCS; + if (ss>=0 && ast_ && ast_[ss]) + if (astWCSIsASkyFrame(ast_[ss])) + return 1; + + return 0; +} + #else + +int FitsImage::hasWCS(Coord::CoordSystem sys) +{ + astClearStatus; + astBegin; + + int ss = sys-Coord::WCS; + if (ss<0 || !newast_) + return 0; + + int rr =0; + int nn = astGetI(newast_,"nframe"); + char cc = ' '; + if (ss) + cc = ss+'@'; + + for (int ii=0; ii<nn; ii++) { + const char* id = astGetC(astGetFrame(newast_,ii+1),"Ident"); + if (cc == id[0]) { + rr =1; + break; + } + } + + astEnd; // now, clean up memory + return rr; +} + int FitsImage::hasWCSEqu(Coord::CoordSystem sys) { astClearStatus; + astBegin; int ss = sys-Coord::WCS; - if (ss>=0 && ast_ && ast_[ss]) - if (astWCSIsASkyFrame(ast_[ss])) { + if (ss<0 || !newast_) + return 0; + + int rr =0; + int nn = astGetI(newast_, "nframe"); + char cc = ' '; + if (ss) + cc = ss+'@'; + + for (int ii=0; ii<nn; ii++) { + AstFrame* ff = (AstFrame*)astGetFrame(newast_,ii+1); + const char* id = astGetC(ff, "Ident"); + if (cc == id[0]) { + + int naxes = astGetI(ff, "Naxes"); + switch (naxes) { + case 2: + rr = astIsASkyFrame(ff); + break; + case 3: + case 4: + { + char* domain = (char*)astGetC(ff, "Domain"); + char* sky = strstr(domain, "SKY"); + rr = sky ? 1 : 0; + } + break; + default: + rr =0; + break; + } + // check for xLON/xLAT and xxLN/xxLT // but GLON/GLAT is ok - const char* str = astGetC(ast_[ss], "System"); + const char* str = astGetC(ff, "System"); if (!strncmp(str,"Unknown",7)) - return 0; - return 1; + rr =0; + break; } + } - return 0; + astEnd; // now, clean up memory + return rr; } -#endif - int FitsImage::hasWCSCel(Coord::CoordSystem sys) { astClearStatus; + astBegin; int ss = sys-Coord::WCS; - if (ss>=0 && ast_ && ast_[ss]) - if (astWCSIsASkyFrame(ast_[ss])) - return 1; + if (ss<0 || !newast_) + return 0; + + int rr =0; + int nn = astGetI(newast_, "nframe"); + char cc = ' '; + if (ss) + cc = ss+'@'; - return 0; + for (int ii=0; ii<nn; ii++) { + AstFrame* ff = (AstFrame*)astGetFrame(newast_,ii+1); + const char* id = astGetC(ff, "Ident"); + if (cc == id[0]) { + + int naxes = astGetI(ff, "Naxes"); + switch (naxes) { + case 2: + rr = astIsASkyFrame(ff); + break; + case 3: + case 4: + { + char* domain = (char*)astGetC(ff, "Domain"); + char* sky = strstr(domain, "SKY"); + rr = sky ? 1 : 0; + } + break; + default: + rr =0; + break; + } + break; + } + } + + astEnd; // now, clean up memory + return rr; } +#endif + // WCSX #ifndef NEWWCS @@ -3458,15 +3583,6 @@ int FitsImage::hasWCS3D(Coord::CoordSystem sys, int aa) { int ss = sys-Coord::WCS; return (aa>=2&&aa<FTY_MAXAXES && sys>=Coord::WCS && wcsx_[ss]) ? 1 : 0; - - /* - if (ss>=0 && ast_ && ast_[ss]) { - int naxes = astGetI(ast_[ss],"Naxes"); - return (aa>=2 && aa<FTY_MAXAXES && naxes >= 3) ? 1 : 0; - } - else - return 0; - */ } double FitsImage::pix2wcsx(double in, Coord::CoordSystem sys, int aa) @@ -3617,6 +3733,63 @@ void FitsImage::astinit(int ss, FitsHead* hd, FitsHead* prim) setAstWCSSkyFrame(ast_[ss],Coord::FK5); } +#ifdef NEWWCS +void FitsImage::astinit(FitsHead* hd, FitsHead* prim) +{ + // just in case + if (!hd) + return; + + newast_ = fits2ast(hd); + if (!newast_) + return; + + astClearStatus; // just to make sure + astBegin; // start memory management + + int naxes = astGetI(newast_,"Naxes"); + switch (naxes) { + case 1: + break; + case 2: + if (astIsASkyFrame(astGetFrame(newast_,AST__CURRENT)) && + astGetI(newast_,"LatAxis") == 1) { + int orr[] = {2,1}; + astPermAxes(newast_,orr); + } + break; + case 3: + case 4: + { + AstFrameSet* ast = newast_; + + int pickc[2] = {1,2}; + AstMapping** mapc = NULL; + AstFrame* permc = (AstFrame*)astPickAxes(ast, 2, pickc, &mapc); + astAddFrame(ast, AST__CURRENT, mapc, permc); + + int isky = astGetI(ast, "Current"); + int pickb[4] = {1, 2, 0, 0}; + AstMapping* mapb; + AstFrame* foo = astFrame(2,"Domain=DATA"); + astPickAxes(foo, naxes, pickb, &mapb); + astInvert(mapb); + astAddFrame(ast, AST__BASE, mapb, foo); + int idata = astGetI(ast, "Current"); + astSetI(ast, "Current", isky); + astSetI(ast, "Base", idata); + } + break; + } + + astEnd; // now, clean up memory + + // set default skyframe + if (astWCSIsASkyFrame(newast_)) + setAstWCSSkyFrame(newast_,Coord::FK5); +} +#endif + void FitsImage::astinit0(int ss, FitsHead* hd, FitsHead* prim) { if (!wcs_[ss]) { diff --git a/tksao/frame/fitsimage.h b/tksao/frame/fitsimage.h index b8fc285..c14c193 100644 --- a/tksao/frame/fitsimage.h +++ b/tksao/frame/fitsimage.h @@ -114,6 +114,9 @@ class FitsImage { WorldCoor** wcs_; // wcs list WCSx** wcsx_; // xth Axis WCS AstFrameSet** ast_; // ast frameset; +#ifdef NEWWCS + AstFrameSet* newast_; // ast frameset; +#endif FitsHead* wcsHeader_; // alt wcs header FitsHead* altHeader_; // wcs header for wfpc2 @@ -141,6 +144,9 @@ class FitsImage { void wcsShow(WorldCoor*); void astinit(int, FitsHead*, FitsHead*); +#ifdef NEWWCS + void astinit(FitsHead*, FitsHead*); +#endif void astinit0(int, FitsHead*, FitsHead*); int checkAstWCS(double, double); AstFrameSet* fits2ast(FitsHead*); @@ -407,9 +413,9 @@ class FitsImage { {return (ast_ && ast_[sys-Coord::WCS]) ? ast_[sys-Coord::WCS] : NULL;} int hasWCS(Coord::CoordSystem); - int hasWCS3D(Coord::CoordSystem, int); int hasWCSEqu(Coord::CoordSystem); int hasWCSCel(Coord::CoordSystem); + int hasWCS3D(Coord::CoordSystem, int); void updateMatrices(Matrix&, Matrix&, Matrix&, Matrix&, Matrix&); void updateMatrices(Matrix3d&, Matrix3d&, Matrix3d&, Matrix3d&); |