From 6d8c1064b17824fa67a081dc3a7efe630677fc58 Mon Sep 17 00:00:00 2001 From: William Joye Date: Wed, 29 Nov 2017 17:16:16 -0500 Subject: update AST WCS --- tksao/frame/fitsimage.C | 91 ++++++++++++++++++++++++++++++++++++++++++++++--- tksao/frame/fitsimage.h | 11 ++++++ tksao/frame/fitsmap.C | 38 ++++++++++++++++++++- 3 files changed, 134 insertions(+), 6 deletions(-) diff --git a/tksao/frame/fitsimage.C b/tksao/frame/fitsimage.C index 49eb5df..ec1387e 100644 --- a/tksao/frame/fitsimage.C +++ b/tksao/frame/fitsimage.C @@ -105,6 +105,9 @@ FitsImage::FitsImage(Context* cx, Tcl_Interp* pp) imageToData3d = Translate3d(-.5, -.5, -.5); dataToImage3d = Translate3d( .5, .5, .5); + imageToRef3d = imageToData3d; + refToImage3d = dataToImage3d; + manageWCS_ =1; wcsx_ =NULL; #ifndef NEWWCS @@ -3119,7 +3122,7 @@ Vector FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys, return Vector(); } #else -Vector FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys, +Vector FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys, Coord::SkyFrame sky) { astClearStatus; // just to make sure @@ -3140,10 +3143,32 @@ Vector FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys, return Vector(); } + +Vector3d FitsImage::pix2wcs(Vector3d in, Coord::CoordSystem sys, + Coord::SkyFrame sky) +{ + astClearStatus; // just to make sure + + if (!hasWCS(sys)) + return Vector(); + + setWCSSystem(newast_,sys); + setWCSSkyFrame(newast_,sky); + + Vector out = wcsTran(newast_, in, 1); + if (astOK && checkWCS(out)) { + if (wcsIsASkyFrame(newast_)) + return out.radToDeg(); + else + return out; + } + + return Vector3d(); +} #endif #ifndef NEWWCS -char* FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys, +char* FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys, Coord::SkyFrame sky, Coord::SkyFormat format, char* lbuf) { @@ -3204,7 +3229,7 @@ char* FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys, return lbuf; } #else -char* FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys, +char* FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys, Coord::SkyFrame sky, Coord::SkyFormat format, char* lbuf) { @@ -3266,7 +3291,7 @@ char* FitsImage::pix2wcs(Vector in, Coord::CoordSystem sys, #endif #ifndef NEWWCS -Vector FitsImage::wcs2pix(Vector in, Coord::CoordSystem sys, +Vector FitsImage::wcs2pix(Vector in, Coord::CoordSystem sys, Coord::SkyFrame sky) { astClearStatus; @@ -3290,7 +3315,7 @@ Vector FitsImage::wcs2pix(Vector in, Coord::CoordSystem sys, return Vector(); } #else -Vector FitsImage::wcs2pix(Vector in, Coord::CoordSystem sys, +Vector FitsImage::wcs2pix(Vector in, Coord::CoordSystem sys, Coord::SkyFrame sky) { astClearStatus; // just to make sure @@ -3310,6 +3335,26 @@ Vector FitsImage::wcs2pix(Vector in, Coord::CoordSystem sys, maperr =1; return Vector(); } + +Vector3d FitsImage::wcs2pix(Vector3d in, Coord::CoordSystem sys, + Coord::SkyFrame sky) +{ + astClearStatus; // just to make sure + + if (hasWCS(sys)) { + setWCSSystem(newast_,sys); + setWCSSkyFrame(newast_,sky); + + if (wcsIsASkyFrame(newast_)) + in *= M_PI/180.; + + Vector out = wcsTran(newast_, in, 0); + if (astOK && checkWCS(out)) + return out; + } + + return Vector3d(); +} #endif #ifndef NEWWCS @@ -4097,6 +4142,42 @@ void FitsImage::wcsTran(AstFrameSet* ast, int npoint, break; } } + +Vector3d FitsImage::wcsTran(AstFrameSet* ast, Vector3d& in, int forward) +{ + int naxes = astGetI(ast,"Naxes"); + switch (naxes) { + case 1: + case 2: + // error + break; + case 3: + { + double pin[3]; + double pout[3]; + pin[0] = in[0]; + pin[1] = in[1]; + pin[2] = in[3]; + astTranN(ast, 1, 3, 1, pin, forward, 3, 1, pout); + return Vector3d(pout[0],pout[1],pout[2]); + } + break; + case 4: + { + double pin[4]; + double pout[4]; + pin[0] = in[0]; + pin[1] = in[1]; + pin[2] = in[3]; + pin[3] = forward ? context_->slice(3) : 0; + astTranN(ast, 1, 4, 1, pin, forward, 4, 1, pout); + return Vector3d(pout[0],pout[1],pout[2]); + } + break; + } + return Vector3d(); +} + #endif #ifndef NEWWCS diff --git a/tksao/frame/fitsimage.h b/tksao/frame/fitsimage.h index 38872f6..44ea871 100644 --- a/tksao/frame/fitsimage.h +++ b/tksao/frame/fitsimage.h @@ -209,6 +209,8 @@ class FitsImage { Matrix3d dataToImage3d; Matrix3d imageToData3d; + Matrix3d refToImage3d; + Matrix3d imageToRef3d; Matrix3d dataToRef3d; Matrix3d refToData3d; @@ -366,6 +368,10 @@ class FitsImage { Vector wcs2pix(Vector, Coord::CoordSystem, Coord::SkyFrame); +#ifdef NEWWCS + Vector3d pix2wcs(Vector3d, Coord::CoordSystem, Coord::SkyFrame); + Vector3d wcs2pix(Vector3d, Coord::CoordSystem, Coord::SkyFrame); +#endif double pix2wcsx(double, Coord::CoordSystem, int); double wcs2pixx(double, Coord::CoordSystem, int); @@ -403,6 +409,7 @@ class FitsImage { void wcsTran(AstFrameSet*, int, Vector*, int, Vector*); double wcsDistance(AstFrameSet*, Vector, Vector); #ifdef NEWWCS + Vector3d wcsTran(AstFrameSet*, Vector3d&, int); double wcsAngle(AstFrameSet*, Vector&, Vector&, Vector&); double wcsAxAngle(AstFrameSet*, Vector&, Vector&); #endif @@ -444,6 +451,10 @@ class FitsImage { Vector mapFromRef(const Vector&, Coord::CoordSystem, Coord::SkyFrame =Coord::FK5); void mapFromRef(const Vector&, Coord::CoordSystem, Coord::SkyFrame, Coord::SkyFormat, char*); Vector mapToRef(const Vector&, Coord::CoordSystem, Coord::SkyFrame =Coord::FK5); +#ifdef NEWWCS + Vector3d mapFromRef(const Vector3d&, Coord::CoordSystem, Coord::SkyFrame =Coord::FK5); + Vector3d mapToRef(const Vector3d&, Coord::CoordSystem, Coord::SkyFrame =Coord::FK5); +#endif double mapFromRef3axis(double, Coord::CoordSystem, int); double mapToRef3axis(double, Coord::CoordSystem, int); double mapLenFromRef(double, Coord::CoordSystem, Coord::DistFormat =Coord::DEGREE); diff --git a/tksao/frame/fitsmap.C b/tksao/frame/fitsmap.C index 60a453a..618aad7 100644 --- a/tksao/frame/fitsmap.C +++ b/tksao/frame/fitsmap.C @@ -20,7 +20,7 @@ Vector FitsImage::mapFromRef(const Vector& vv, Coord::CoordSystem out, return vv * refToDetector; default: if (hasWCS(out)) - return pix2wcs(vv * refToImage, out, sky); + return pix2wcs(vv * refToImage3d, out, sky); } return Vector(); @@ -60,6 +60,42 @@ Vector FitsImage::mapToRef(const Vector& vv, Coord::CoordSystem in, return Vector(); } +#ifdef NEWWCS +Vector3d FitsImage::mapFromRef(const Vector3d& vv, Coord::CoordSystem out, + Coord::SkyFrame sky) +{ + switch (out) { + case Coord::IMAGE: + case Coord::PHYSICAL: + case Coord::AMPLIFIER: + case Coord::DETECTOR: + return vv * refToImage3d; + default: + if (hasWCS(out)) + return pix2wcs(vv * refToImage3d, out, sky); + } + + return Vector3d(); +} + +Vector3d FitsImage::mapToRef(const Vector3d& vv, Coord::CoordSystem in, + Coord::SkyFrame sky) +{ + switch (in) { + case Coord::IMAGE: + case Coord::PHYSICAL: + case Coord::AMPLIFIER: + case Coord::DETECTOR: + return vv * imageToRef3d; + default: + if (hasWCS(in)) + return wcs2pix(vv, in, sky) * imageToRef3d; + } + + return Vector3d(); +} +#endif + void FitsImage::listFromRef(ostream& str, const Vector& vv, Coord::CoordSystem sys, Coord::SkyFrame sky, Coord::SkyFormat format) -- cgit v0.12