From 5e9b73e09b0483a8b92215b1ec2d0893de730232 Mon Sep 17 00:00:00 2001
From: William Joye <wjoye@cfa.harvard.edu>
Date: Thu, 9 Nov 2017 14:30:29 -0500
Subject: update AST WCS
---
tksao/frame/fitsimage.C | 60 ++++++++++++++++++++++++++++++++++++++++---------
tksao/frame/fitsimage.h | 6 ++++-
tksao/frame/grid2d.C | 55 ++++++++++++++++++++++++++++++++++++++++-----
3 files changed, 105 insertions(+), 16 deletions(-)
diff --git a/tksao/frame/fitsimage.C b/tksao/frame/fitsimage.C
index 81ff0a7..0c7eba5 100644
--- a/tksao/frame/fitsimage.C
+++ b/tksao/frame/fitsimage.C
@@ -3346,7 +3346,9 @@ char* FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys,
return lbuf;
}
+
#else
+
char* FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys,
Coord::SkyFrame sky, Coord::SkyFormat format,
char* lbuf)
@@ -3925,8 +3927,8 @@ void FitsImage::astinit(int ss, FitsHead* hd, FitsHead* prim)
setAstWCSSystem(ast_[ss], (Coord::CoordSystem)(ss+Coord::WCS));
- astClearStatus; // just to make sure
- astBegin; // start memory management
+ // astClearStatus; // just to make sure
+ // astBegin; // start memory management
int naxes = astGetI(ast_[ss],"Naxes");
switch (naxes) {
@@ -3942,6 +3944,7 @@ void FitsImage::astinit(int ss, FitsHead* hd, FitsHead* prim)
case 3:
case 4:
{
+ if (0) {
AstFrameSet* ast = ast_[ss];
int pickc[2] = {1,2};
@@ -3959,11 +3962,12 @@ void FitsImage::astinit(int ss, FitsHead* hd, FitsHead* prim)
int idata = astGetI(ast, "Current");
astSetI(ast, "Current", isky);
astSetI(ast, "Base", idata);
+ }
}
break;
}
- astEnd; // now, clean up memory
+ // astEnd; // now, clean up memory
#endif
// set default skyframe
@@ -3982,9 +3986,6 @@ void FitsImage::astinit(FitsHead* hd, FitsHead* prim)
if (!newast_)
return;
- astClearStatus; // just to make sure
- astBegin; // start memory management
-
int naxes = astGetI(newast_,"Naxes");
switch (naxes) {
case 1:
@@ -3999,6 +4000,10 @@ void FitsImage::astinit(FitsHead* hd, FitsHead* prim)
case 3:
case 4:
{
+ if (0) {
+ astClearStatus; // just to make sure
+ astBegin; // start memory management
+
AstFrameSet* ast = newast_;
int pickc[2] = {1,2};
@@ -4016,13 +4021,15 @@ void FitsImage::astinit(FitsHead* hd, FitsHead* prim)
int idata = astGetI(ast, "Current");
astSetI(ast, "Current", isky);
astSetI(ast, "Base", idata);
+
+ astEnd; // now, clean up memory
+ }
}
break;
}
- astEnd; // now, clean up memory
-
- setAstWCSSkyFrame(newast_,Coord::FK5);
+ if (astWCSIsASkyFrame(newast_))
+ setAstWCSSkyFrame(newast_,Coord::FK5);
}
#endif
@@ -4205,7 +4212,40 @@ void FitsImage::astWCSTran(AstFrameSet* ast, int npoint,
int forward,
double* xout, double* yout)
{
- astTran2(ast, npoint, xin, yin, forward, xout, yout);
+ // astTran2(ast, npoint, xin, yin, forward, xout, yout);
+
+ int naxes = astGetI(ast,"Naxes");
+ switch (naxes) {
+ case 1:
+ // error
+ break;
+ case 2:
+ astTran2(ast, npoint, xin, yin, forward, xout, yout);
+ break;
+ case 3:
+ {
+ double* ptr_in[3];
+ ptr_in[0] = (double*)xin;
+ ptr_in[1] = (double*)yin;
+ ptr_in[2] = new double[npoint];
+ for (int kk=0; kk<npoint; kk++)
+ ptr_in[2][kk] = 1;
+
+
+ double* ptr_out[3];
+ ptr_out[0] = (double*)xout;
+ ptr_out[1] = (double*)yout;
+ ptr_out[2] = new double[npoint];
+
+ astTranP(ast, npoint, 3, (const double**)ptr_in, forward, 3, ptr_out);
+
+ if (ptr_in[2])
+ delete [] ptr_in[2];
+ if (ptr_out[2])
+ delete [] ptr_out[2];
+ }
+ break;
+ }
}
#endif
diff --git a/tksao/frame/fitsimage.h b/tksao/frame/fitsimage.h
index c14c193..3759285 100644
--- a/tksao/frame/fitsimage.h
+++ b/tksao/frame/fitsimage.h
@@ -409,9 +409,13 @@ class FitsImage {
#endif
void setAstWCSSkyFrame(AstFrameSet*, Coord::SkyFrame);
void setAstWCSFormat(AstFrameSet*, int, const char*);
+#ifndef NEWWCS
AstFrameSet* getAST(Coord::CoordSystem sys)
{return (ast_ && ast_[sys-Coord::WCS]) ? ast_[sys-Coord::WCS] : NULL;}
-
+#else
+ AstFrameSet* getAST(Coord::CoordSystem sys) {return newast_;}
+#endif
+
int hasWCS(Coord::CoordSystem);
int hasWCSEqu(Coord::CoordSystem);
int hasWCSCel(Coord::CoordSystem);
diff --git a/tksao/frame/grid2d.C b/tksao/frame/grid2d.C
index 6d5bf22..aea9320 100644
--- a/tksao/frame/grid2d.C
+++ b/tksao/frame/grid2d.C
@@ -68,9 +68,9 @@ int Grid2d::doit(RenderMode rm)
break;
default:
{
+#ifndef NEWWCS
AstFrameSet* ast = (AstFrameSet*)astCopy(fits->getAST(system_));
-#ifndef NEWWCS
// set desired skyformat
if (fits->astWCSIsASkyFrame(ast))
fits->setAstWCSSkyFrame(ast, sky_);
@@ -81,14 +81,59 @@ int Grid2d::doit(RenderMode rm)
astInvert(ast);
astAddFrame(frameSet,2,astUnitMap(2,""),ast);
astSetI(frameSet,"current",astGetI(frameSet,"nframe"));
-
#else
+ astClearStatus; // just to make sure
+ astBegin; // start memory management
+
+ AstFrameSet* ast = (AstFrameSet*)astCopy(fits->getAST(system_));
+
+ fits->setAstWCSSystem(ast, system_);
+ fits->setAstWCSSkyFrame(ast, sky_);
+
+ int offset =0;
+ int naxes = astGetI(ast,"Naxes");
+ switch (naxes) {
+ case 1:
+ case 2:
+ break;
+ case 3:
+ case 4:
+ {
+ int pick[2] = {1, 2};
+ astInvert(ast);
+ AstMapping* mapb =NULL;
+ AstFrame* permb = (AstFrame*)astPickAxes(ast, 2, pick, &mapb);
+ astAddFrame(ast, AST__CURRENT, mapb, permb);
+ astInvert(ast);
+
+ AstMapping** mapc =NULL;
+ AstFrame* permc = (AstFrame*)astPickAxes(ast, 2, pick, &mapc);
+ astAddFrame(ast, AST__CURRENT, mapc, permc);
+
+ if (0) {
+ 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);
+ offset =1;
+ }
+ }
+ break;
+ }
+
// add wcs to frameset
// this will link frameset to wcs with unitMap
astInvert(ast);
- astAddFrame(frameSet,AST__CURRENT,astUnitMap(2,""),ast);
- fits->setAstWCSSystem(frameSet,system_);
- fits->setAstWCSSkyFrame(frameSet,sky_);
+ astAddFrame(frameSet, AST__CURRENT, astUnitMap(2,""), ast);
+ astSetI(frameSet,"Current",astGetI(frameSet,"nframe")-offset);
+
+ astEnd; // now, clean up memory
#endif
}
}
--
cgit v0.12