From 4294afb4810615addc9ab815489f008a2dcce97b Mon Sep 17 00:00:00 2001 From: William Joye Date: Sun, 3 Dec 2017 15:57:22 -0500 Subject: update AST WCS --- tksao/frame/base.C | 2 +- tksao/frame/basecommand.C | 18 +++---- tksao/frame/fitsimage.C | 120 ++++++++++++++++++++++------------------------ tksao/frame/fitsimage.h | 16 +++---- tksao/frame/fitsmap.C | 8 ++-- tksao/frame/frame3dbase.C | 2 +- tksao/frame/frblt.C | 2 +- tksao/frame/grid3d.C | 10 ++-- 8 files changed, 85 insertions(+), 93 deletions(-) diff --git a/tksao/frame/base.C b/tksao/frame/base.C index ebf8a9b..594573a 100644 --- a/tksao/frame/base.C +++ b/tksao/frame/base.C @@ -916,7 +916,7 @@ int Base::hasWCSCel(Coord::CoordSystem sys) int Base::hasWCS3D(Coord::CoordSystem sys) { - return currentContext->cfits && currentContext->cfits->hasWCS3D(sys,2); + return currentContext->cfits && currentContext->cfits->hasWCS3D(sys); } Vector Base::imageCenter(FrScale::SecMode mode) diff --git a/tksao/frame/basecommand.C b/tksao/frame/basecommand.C index e466e11..e39c90d 100644 --- a/tksao/frame/basecommand.C +++ b/tksao/frame/basecommand.C @@ -568,8 +568,8 @@ void Base::crop3dCmd(double z0, double z1, Coord::CoordSystem sys) return; // ff/tt ranges 0-n - double ff = ptr->mapToRef3axis(z0,sys,2); - double tt = ptr->mapToRef3axis(z1,sys,2); + double ff = ptr->mapToRef3axis(z0,sys); + double tt = ptr->mapToRef3axis(z1,sys); // params is a BBOX in DATA coords 0-n currentContext->setCrop3dParams(ff-.5,tt+.5); @@ -1384,8 +1384,8 @@ void Base::getCoord3axisCmd(double vv, Coord::CoordSystem in, printDouble(vv); else { // use first slice - double rr = currentContext->fits->mapToRef3axis(vv,in,ss); - double tt = currentContext->fits->mapFromRef3axis(rr,out,ss); + double rr = currentContext->fits->mapToRef3axis(vv,in); + double tt = currentContext->fits->mapFromRef3axis(rr,out); printDouble(tt); } } @@ -1440,8 +1440,8 @@ void Base::getCrop3dCmd(Coord::CoordSystem sys) FitsZBound* zparams = currentContext->getDataParams(currentContext->secMode()); - double ff = ptr->mapFromRef3axis(zparams->zmin+.5,sys,2); - double tt = ptr->mapFromRef3axis(zparams->zmax-.5,sys,2); + double ff = ptr->mapFromRef3axis(zparams->zmin+.5,sys); + double tt = ptr->mapFromRef3axis(zparams->zmax-.5,sys); ostringstream str; str << ff << ' ' << tt << ends; @@ -1793,7 +1793,7 @@ void Base::getFitsSliceCmd(int id, Coord::CoordSystem sys) { if (currentContext->fits) { int ss = currentContext->slice(id); - double rr = currentContext->fits->mapFromRef3axis(ss,sys,id); + double rr = currentContext->fits->mapFromRef3axis(ss,sys); printDouble(rr); } else @@ -2826,8 +2826,8 @@ void Base::sliceCmd(int id, int ss) void Base::sliceCmd(int id, double vv, Coord::CoordSystem sys) { - double rr = currentContext->fits->mapToRef3axis(vv,sys,id); - int ss = currentContext->fits->mapFromRef3axis(rr,Coord::IMAGE,id); + double rr = currentContext->fits->mapToRef3axis(vv,sys); + int ss = currentContext->fits->mapFromRef3axis(rr,Coord::IMAGE); // IMAGE (ranges 1-n) setSlice(id,ss); diff --git a/tksao/frame/fitsimage.C b/tksao/frame/fitsimage.C index f9fe867..04713e8 100644 --- a/tksao/frame/fitsimage.C +++ b/tksao/frame/fitsimage.C @@ -53,11 +53,9 @@ extern "C" { WCSx::WCSx() { - for (int ii=0; iifind(scrpix) && hd->find(scrval)) { - if (!wcsx_[ii]) - wcsx_[ii] = new WCSx(); - wcsx_[ii]->crpix[jj] = hd->getReal(scrpix,0); - wcsx_[ii]->crval[jj] = hd->getReal(scrval,0); - - float cd = hd->getReal(scd,0); - float pc = hd->getReal(spc,0); - float cdelt = hd->getReal(scdelt,0); - if (cd) - wcsx_[ii]->cd[jj] = cd; - else if (pc && cdelt) - wcsx_[ii]->cd[jj] = pc * cdelt; - else if (cdelt) - wcsx_[ii]->cd[jj] = cdelt; - else - wcsx_[ii]->cd[jj] = 1; - } + scrpix[5] = '1'+jj; + scrpix[6] = !ii ? ' ' : '@'+jj; + scrval[5] = '1'+jj; + scrval[6] = !ii ? ' ' : '@'+jj; + scd[2] = '1'+jj; + scd[4] = '1'+jj; + scd[5] = !ii ? ' ' : '@'+jj; + spc[2] = '1'+jj; + spc[4] = '1'+jj; + spc[5] = !ii ? ' ' : '@'+jj; + scdelt[5] = '1'+jj; + scdelt[6] = !ii ? ' ' : '@'+jj; + + if (hd->find(scrpix) && hd->find(scrval)) { + if (!wcsx_[ii]) + wcsx_[ii] = new WCSx(); + wcsx_[ii]->crpix = hd->getReal(scrpix,0); + wcsx_[ii]->crval = hd->getReal(scrval,0); + + float cd = hd->getReal(scd,0); + float pc = hd->getReal(spc,0); + float cdelt = hd->getReal(scdelt,0); + if (cd) + wcsx_[ii]->cd = cd; + else if (pc && cdelt) + wcsx_[ii]->cd = pc * cdelt; + else if (cdelt) + wcsx_[ii]->cd = cdelt; + else + wcsx_[ii]->cd = 1; } } @@ -1273,18 +1270,13 @@ void FitsImage::initWCS() if (DebugWCS) { for (int ii=0; iicd[jj]) { - cerr << "wcsx" << (char)(!ii ? ' ' : '@'+ii) - << "[" << ii << "]->crpix[" << jj << "]=" - << wcsx_[ii]->crpix[jj] << endl; - cerr << "wcsx" << (char)(!ii ? ' ' : '@'+ii) - << "[" << ii << "]->crval[" << jj << "]=" - << wcsx_[ii]->crval[jj] << endl; - cerr << "wcsx" << (char)(!ii ? ' ' : '@'+ii) - << "[" << ii << "]->cd[" << jj << "]=" - << wcsx_[ii]->cd[jj] << endl; - } + if (wcsx_[ii]->cd) { + cerr << "wcsx" << (char)(!ii ? ' ' : '@'+ii) + << "[" << ii << "]->crpix=" << wcsx_[ii]->crpix << endl; + cerr << "wcsx" << (char)(!ii ? ' ' : '@'+ii) + << "[" << ii << "]->crval=" << wcsx_[ii]->crval << endl; + cerr << "wcsx" << (char)(!ii ? ' ' : '@'+ii) + << "[" << ii << "]->cd=" << wcsx_[ii]->cd << endl; } } } @@ -3537,27 +3529,27 @@ int FitsImage::hasWCSCel(Coord::CoordSystem sys) #ifndef NEWWCS -int FitsImage::hasWCS3D(Coord::CoordSystem sys, int aa) +int FitsImage::hasWCS3D(Coord::CoordSystem sys) { int ss = sys-Coord::WCS; - return (aa>=2&&aa=Coord::WCS && wcsx_[ss]) ? 1 : 0; + return (sys>=Coord::WCS && wcsx_[ss]) ? 1 : 0; } -double FitsImage::pix2wcsx(double in, Coord::CoordSystem sys, int aa) +double FitsImage::pix2wcsx(double in, Coord::CoordSystem sys) { - if (hasWCS3D(sys,aa)) { + if (hasWCS3D(sys)) { int ss = sys-Coord::WCS; - return (in-wcsx_[ss]->crpix[aa])*wcsx_[ss]->cd[aa] + wcsx_[ss]->crval[aa]; + return (in-wcsx_[ss]->crpix)*wcsx_[ss]->cd + wcsx_[ss]->crval; } else return in; } -double FitsImage::wcs2pixx(double in, Coord::CoordSystem sys, int aa) +double FitsImage::wcs2pixx(double in, Coord::CoordSystem sys) { - if (hasWCS3D(sys,aa)) { + if (hasWCS3D(sys)) { int ss = sys-Coord::WCS; - return (in-wcsx_[ss]->crval[aa])/wcsx_[ss]->cd[aa] + wcsx_[ss]->crpix[aa]; + return (in-wcsx_[ss]->crval)/wcsx_[ss]->cd + wcsx_[ss]->crpix; } else return in; @@ -3565,7 +3557,7 @@ double FitsImage::wcs2pixx(double in, Coord::CoordSystem sys, int aa) #else -int FitsImage::hasWCS3D(Coord::CoordSystem sys, int aa) +int FitsImage::hasWCS3D(Coord::CoordSystem sys) { if (!newast_ || syscrpix[aa])*wcsx_[ss]->cd[aa] + wcsx_[ss]->crval[aa]; + return (in-wcsx_[ss]->crpix)*wcsx_[ss]->cd + wcsx_[ss]->crval; } else return in; } -double FitsImage::wcs2pixx(double in, Coord::CoordSystem sys, int aa) +double FitsImage::wcs2pixx(double in, Coord::CoordSystem sys) { - if (hasWCS3D(sys,aa)) { + if (hasWCS3D(sys)) { int ss = sys-Coord::WCS; - return (in-wcsx_[ss]->crval[aa])/wcsx_[ss]->cd[aa] + wcsx_[ss]->crpix[aa]; + return (in-wcsx_[ss]->crval)/wcsx_[ss]->cd + wcsx_[ss]->crpix; } else return in; diff --git a/tksao/frame/fitsimage.h b/tksao/frame/fitsimage.h index 09de2f3..b7713bd 100644 --- a/tksao/frame/fitsimage.h +++ b/tksao/frame/fitsimage.h @@ -45,9 +45,9 @@ extern "C" { class WCSx { public: - double crpix[FTY_MAXAXES]; - double crval[FTY_MAXAXES]; - double cd[FTY_MAXAXES]; + double crpix; + double crval; + double cd; public: WCSx(); @@ -373,8 +373,8 @@ class FitsImage { Vector3d pix2wcs(const Vector3d&, Coord::CoordSystem, Coord::SkyFrame); Vector3d wcs2pix(const Vector3d&, Coord::CoordSystem, Coord::SkyFrame); #endif - double pix2wcsx(double, Coord::CoordSystem, int); - double wcs2pixx(double, Coord::CoordSystem, int); + double pix2wcsx(double, Coord::CoordSystem); + double wcs2pixx(double, Coord::CoordSystem); void altWCS(istream&); void appendWCS(istream&); @@ -430,7 +430,7 @@ class FitsImage { int hasWCS(Coord::CoordSystem); int hasWCSEqu(Coord::CoordSystem); int hasWCSCel(Coord::CoordSystem); - int hasWCS3D(Coord::CoordSystem, int); + int hasWCS3D(Coord::CoordSystem); void updateMatrices(Matrix&, Matrix&, Matrix&, Matrix&, Matrix&); void updateMatrices(Matrix3d&, Matrix3d&, Matrix3d&, Matrix3d&); @@ -456,8 +456,8 @@ class FitsImage { 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 mapFromRef3axis(double, Coord::CoordSystem); + double mapToRef3axis(double, Coord::CoordSystem); double mapLenFromRef(double, Coord::CoordSystem, Coord::DistFormat =Coord::DEGREE); Vector mapLenFromRef(const Vector&, Coord::CoordSystem, Coord::DistFormat =Coord::DEGREE); double mapLenToRef(double, Coord::CoordSystem, Coord::DistFormat =Coord::DEGREE); diff --git a/tksao/frame/fitsmap.C b/tksao/frame/fitsmap.C index 2e28382..17d2369 100644 --- a/tksao/frame/fitsmap.C +++ b/tksao/frame/fitsmap.C @@ -518,7 +518,7 @@ void FitsImage::listDistFromRef(ostream& str, // 3D -double FitsImage::mapFromRef3axis(double vv, Coord::CoordSystem out, int ss) +double FitsImage::mapFromRef3axis(double vv, Coord::CoordSystem out) { switch (out) { case Coord::IMAGE: @@ -527,11 +527,11 @@ double FitsImage::mapFromRef3axis(double vv, Coord::CoordSystem out, int ss) case Coord::DETECTOR: return vv+.5; default: - return pix2wcsx(vv+.5,out,ss); + return pix2wcsx(vv+.5,out); } } -double FitsImage::mapToRef3axis(double vv, Coord::CoordSystem in, int ss) +double FitsImage::mapToRef3axis(double vv, Coord::CoordSystem in) { switch (in) { case Coord::IMAGE: @@ -540,6 +540,6 @@ double FitsImage::mapToRef3axis(double vv, Coord::CoordSystem in, int ss) case Coord::DETECTOR: return vv-.5; default: - return wcs2pixx(vv,in,ss) -.5; + return wcs2pixx(vv,in) -.5; } } diff --git a/tksao/frame/frame3dbase.C b/tksao/frame/frame3dbase.C index 0b56286..3b39d3e 100644 --- a/tksao/frame/frame3dbase.C +++ b/tksao/frame/frame3dbase.C @@ -210,7 +210,7 @@ void Frame3dBase::coordToTclArray(FitsImage* ptr, const Vector3d& vv, doubleToTclArray(rr[0], var, base, "x"); doubleToTclArray(rr[1], var, base, "y"); - double ss = ptr->mapFromRef3axis(((Vector3d&)vv)[2],out,2); + double ss = ptr->mapFromRef3axis(((Vector3d&)vv)[2],out); doubleToTclArray(ss, var, base, "z"); } #else diff --git a/tksao/frame/frblt.C b/tksao/frame/frblt.C index 95d812b..1584197 100644 --- a/tksao/frame/frblt.C +++ b/tksao/frame/frblt.C @@ -283,7 +283,7 @@ int Base::markerAnalysisPlot3d(Marker* pp, double** x, double** y, // main loop SETSIGBUS for (int kk=0; kkmapFromRef3axis(kk+.5+zparams->zmin, sys, 2); + (*x)[kk] = ptr->mapFromRef3axis(kk+.5+zparams->zmin, sys); bool* mptr=msk; long* iptr=idx; diff --git a/tksao/frame/grid3d.C b/tksao/frame/grid3d.C index cf090f5..6850c2f 100644 --- a/tksao/frame/grid3d.C +++ b/tksao/frame/grid3d.C @@ -20,13 +20,13 @@ void bar(AstMapping* that, int npoint, int ncoord_in, const double* ptr_in[], if (forward) { for (int ii=0; iicrpix[2]) * - wcsx[0]->cd[2] + wcsx[0]->crval[2]; + (*ptr_out)[ii] = (.5 + (*ptr_in)[ii] - wcsx[0]->crpix) * + wcsx[0]->cd + wcsx[0]->crval; } else { for (int ii=0; iicrval[2]) / - wcsx[0]->cd[2] + wcsx[0]->crpix[2] -.5; + (*ptr_out)[ii] = ((*ptr_in)[ii] - wcsx[0]->crval) / + wcsx[0]->cd + wcsx[0]->crpix -.5; } } @@ -124,7 +124,7 @@ int Grid3d::doit(RenderMode rm) AstFrame* zbase = astFrame(1,""); AstFrame* zcurr = astFrame(1,""); AstMapping* zmap; - if (fits->hasWCS3D(system_,2)) { + if (fits->hasWCS3D(system_)) { astIntraReg("foo",1,1,bar,0,"testing","me","you"); if (!(zmap = (AstMapping*)astIntraMap("foo",1,1,""))) return 0; -- cgit v0.12