From cba9ac56699485e90abad7bdcfa7473600010937 Mon Sep 17 00:00:00 2001 From: William Joye Date: Tue, 6 Nov 2018 16:39:59 -0500 Subject: fix savefits wcs for dss --- tksao/frame/fitsimage.C | 16 ++++ tksao/frame/fitsimage.h | 1 + tksao/frame/frsave.C | 225 ++++++++++++++---------------------------------- 3 files changed, 83 insertions(+), 159 deletions(-) diff --git a/tksao/frame/fitsimage.C b/tksao/frame/fitsimage.C index 889befd..474326d 100644 --- a/tksao/frame/fitsimage.C +++ b/tksao/frame/fitsimage.C @@ -84,6 +84,7 @@ FitsImage::FitsImage(Context* cx, Tcl_Interp* pp) manageWCS_ =1; ast_ =NULL; + encoding_ =NULL; wcs_ =NULL; wcsNaxes_ =NULL; @@ -159,6 +160,8 @@ FitsImage::~FitsImage() if (manageWCS_) { if (ast_) astAnnul(ast_); + if (encoding_) + delete [] encoding_; if (wcs_) delete [] wcs_; @@ -1062,6 +1065,9 @@ void FitsImage::initWCS(FitsHead* hd) if (ast_) astAnnul(ast_); ast_ =NULL; + if (encoding_) + delete [] encoding_; + encoding_ =NULL; if (wcs_) delete [] wcs_; @@ -1104,6 +1110,7 @@ void FitsImage::initWCS(FitsHead* hd) while (sptr) { if (sptr == this) { ast_ = ptr->ast_; + encoding_ = ptr->encoding_; wcs_ = ptr->wcs_; wcsNaxes_ = ptr->wcsNaxes_; @@ -1134,6 +1141,9 @@ void FitsImage::initWCS(FitsHead* hd) if (ast_) astAnnul(ast_); ast_ =NULL; + if (encoding_) + delete [] encoding_; + encoding_ =NULL; ast_ = fits2ast(hd); if (!ast_) @@ -3464,6 +3474,12 @@ AstFrameSet* FitsImage::fits2ast(FitsHead* hd) astClearStatus; } + // encoding + // must come after astPutsFits and before astRead + const char* encode = astGetC(chan, "Encoding"); + if (encode) + encoding_ = dupstr(encode); + // we may have an error, just reset astClearStatus; astClear(chan, "Card"); diff --git a/tksao/frame/fitsimage.h b/tksao/frame/fitsimage.h index 1acc5bf..5c4fe87 100644 --- a/tksao/frame/fitsimage.h +++ b/tksao/frame/fitsimage.h @@ -134,6 +134,7 @@ class FitsImage { public: AstFrameSet* ast_; // ast frameset; + const char* encoding_; // ast encoding of original header private: char* root(const char*); diff --git a/tksao/frame/frsave.C b/tksao/frame/frsave.C index b0f854e..4a3aee1 100644 --- a/tksao/frame/frsave.C +++ b/tksao/frame/frsave.C @@ -376,57 +376,42 @@ void FrameBase::saveFitsResampleKeyword(OutFitsStream& str, FitsHead& dst) FitsHead* src = currentContext->fits->head(); Vector center = Vector(options->width, options->height)/2.; + // center mx + Matrix cc = + Translate(-center) * + Translate(1,0) * + FlipY() * + Translate(center); + // OBJECT char* object = src->getString("OBJECT"); if (object) dst.appendString("OBJECT", object, NULL); - // DATE-OBS - char* date = src->getString("DATE"); - if (date) - dst.appendString("DATE", date, NULL); - - char* dateobs = src->getString("DATE-OBS"); - if (dateobs) - dst.appendString("DATE-OBS", dateobs, NULL); - - char* timeobs = src->getString("TIME-OBS"); - if (timeobs) - dst.appendString("TIME-OBS", timeobs, NULL); - - char* dateend = src->getString("DATE-END"); - if (dateend) - dst.appendString("DATE-END", dateend, NULL); - - char* timeend = src->getString("TIME-END"); - if (timeend) - dst.appendString("TIME-END", timeend, NULL); + // LTMV + if (currentContext->fits->hasLTMV()) { + Matrix ltmv = currentContext->fits->physicalToRef * refToWidget * cc; + dst.appendReal("LTM1_1", ltmv[0][0], 9, NULL); + dst.appendReal("LTM1_2", ltmv[0][1], 9, NULL); + dst.appendReal("LTM2_1", ltmv[1][0], 9, NULL); + dst.appendReal("LTM2_2", ltmv[1][1], 9, NULL); + dst.appendReal("LTV1", ltmv[2][0], 9, NULL); + dst.appendReal("LTV2", ltmv[2][1], 9, NULL); + } - // LTMV,DTMV if (!isMosaic()) { - if (currentContext->fits->hasLTMV()) { - Matrix ltmv = currentContext->fits->physicalToRef * refToWidget * - Translate(-center) * - Translate(1,0) * - FlipY() * - Translate(center); - - dst.appendReal("LTM1_1", ltmv[0][0], 9, NULL); - dst.appendReal("LTM1_2", ltmv[0][1], 9, NULL); - dst.appendReal("LTM2_1", ltmv[1][0], 9, NULL); - dst.appendReal("LTM2_2", ltmv[1][1], 9, NULL); - dst.appendReal("LTV1", ltmv[2][0], 9, NULL); - dst.appendReal("LTV2", ltmv[2][1], 9, NULL); + if (currentContext->fits->hasATMV()) { + Matrix dtmv = currentContext->fits->amplifierToRef * refToWidget * cc; + dst.appendReal("ATM1_1", dtmv[0][0], 9, NULL); + dst.appendReal("ATM1_2", dtmv[0][1], 9, NULL); + dst.appendReal("ATM2_1", dtmv[1][0], 9, NULL); + dst.appendReal("ATM2_2", dtmv[1][1], 9, NULL); + dst.appendReal("ATV1", dtmv[2][0], 9, NULL); + dst.appendReal("ATV2", dtmv[2][1], 9, NULL); } - } - else { - if (currentContext->fits->hasDTMV()) { - Matrix dtmv = currentContext->fits->detectorToRef * refToWidget * - Translate(-center) * - Translate(1,0) * - FlipY() * - Translate(center); + if (currentContext->fits->hasDTMV()) { + Matrix dtmv = currentContext->fits->detectorToRef * refToWidget * cc; dst.appendReal("DTM1_1", dtmv[0][0], 9, NULL); dst.appendReal("DTM1_2", dtmv[0][1], 9, NULL); dst.appendReal("DTM2_1", dtmv[1][0], 9, NULL); @@ -437,124 +422,46 @@ void FrameBase::saveFitsResampleKeyword(OutFitsStream& str, FitsHead& dst) } // WCS - if (currentContext->fits->hasWCS(Coord::WCS)) { - if (src->find("RADESYS")) - dst.appendString("RADESYS", src->getString("RADESYS"), NULL); - if (src->find("EQUINOX")) - dst.appendReal("EQUINOX", src->getReal("EQUINOX",2000), 9, NULL); - if (src->find("EPOCH")) - dst.appendReal("EPOCH", src->getReal("EPOCH",2000), 9, NULL); - if (src->find("MJD-OBS")) - dst.appendReal("MJD-OBS", src->getReal("MJD-OBS",51544), 9, NULL); - if (src->find("CTYPE1")) - dst.appendString("CTYPE1", src->getString("CTYPE1"), NULL); - if (src->find("CTYPE2")) - dst.appendString("CTYPE2", src->getString("CTYPE2"), NULL); - if (src->find("CRVAL1")) - dst.appendReal("CRVAL1", src->getReal("CRVAL1",0), 9, NULL); - if (src->find("CRVAL2")) - dst.appendReal("CRVAL2", src->getReal("CRVAL2",0), 9, NULL); - if (src->find("CUNIT1")) - dst.appendString("CUNIT1", src->getString("CUNIT1"), NULL); - if (src->find("CUNIT2")) - dst.appendString("CUNIT2", src->getString("CUNIT2"), NULL); - - // crpix - if (src->find("CRPIX1") || src->find("CRPIX2")) { - double crpix1 = src->getReal("CRPIX1",0); - double crpix2 = src->getReal("CRPIX2",0); - - Vector crpix = Vector(crpix1,crpix2) * - currentContext->fits->imageToWidget * - Translate(-center) * - Translate(1,0) * - FlipY() * - Translate(center); - - dst.appendReal("CRPIX1", crpix[0], 9, NULL); - dst.appendReal("CRPIX2", crpix[1], 9, NULL); - } - - // cd - if (src->find("CD1_1") || src->find("CD1_2") || - src->find("CD2_1") || src->find("CD2_2")) { - // cd keywords - double cd11 = src->getReal("CD1_1",0); - double cd12 = src->getReal("CD1_2",0); - double cd21 = src->getReal("CD2_1",0); - double cd22 = src->getReal("CD2_2",0); - - Matrix cd = Matrix(cd11, cd12, cd21, cd22,0,0) * - currentContext->fits->imageToRef * refToUser * - wcsOrientationMatrix * - Rotate(wcsRotation) * - orientationMatrix * - Scale(zoom_.invert()) * - Rotate(rotation) * - Translate(center) * - Translate(-center) * - Translate(1,0) * - FlipY() * - Translate(center); - - dst.appendReal("CD1_1", cd.matrix(0,0), 9, NULL); - dst.appendReal("CD1_2", cd.matrix(0,1), 9, NULL); - dst.appendReal("CD2_1", cd.matrix(1,0), 9, NULL); - dst.appendReal("CD2_2", cd.matrix(1,1), 9, NULL); - } - else if (src->find("PC1_1") || src->find("PC1_2") || - src->find("PC2_1") || src->find("PC2_2")) { - // pc keywords - double pc11 = src->getReal("PC1_1",1); - double pc12 = src->getReal("PC1_2",0); - double pc21 = src->getReal("PC2_1",0); - double pc22 = src->getReal("PC2_2",1); - double cdelt1 = src->getReal("CDELT1",1); - double cdelt2 = src->getReal("CDELT2",1); - - Matrix cd = Scale(cdelt1,cdelt2) * Matrix(pc11,pc12,pc21,pc22,0,0) * - currentContext->fits->imageToRef * refToUser * - wcsOrientationMatrix * - Rotate(wcsRotation) * - orientationMatrix * - Scale(zoom_.invert()) * - Rotate(rotation) * - Translate(center) * - Translate(-center) * - Translate(1,0) * - FlipY() * - Translate(center); - - dst.appendReal("CD1_1", cd.matrix(0,0), 9, NULL); - dst.appendReal("CD1_2", cd.matrix(0,1), 9, NULL); - dst.appendReal("CD2_1", cd.matrix(1,0), 9, NULL); - dst.appendReal("CD2_2", cd.matrix(1,1), 9, NULL); - } - else if (src->find("CDELT1") || src->find("CDELT2")) { - // crota2 - double cdelt1 = src->getReal("CDELT1",1); - double cdelt2 = src->getReal("CDELT2",1); - double crot2 = src->getReal("CROT2",0); - - Matrix cd = Scale(cdelt1,cdelt2) * Rotate(crot2) * - currentContext->fits->imageToRef * refToUser * - wcsOrientationMatrix * - Rotate(wcsRotation) * - orientationMatrix * - Scale(zoom_.invert()) * - Rotate(rotation) * - Translate(center) * - Translate(-center) * - Translate(1,0) * - FlipY() * - Translate(center); - - dst.appendReal("CD1_1", cd.matrix(0,0), 9, NULL); - dst.appendReal("CD1_2", cd.matrix(0,1), 9, NULL); - dst.appendReal("CD2_1", cd.matrix(1,0), 9, NULL); - dst.appendReal("CD2_2", cd.matrix(1,1), 9, NULL); - } + astClearStatus; // just to make sure + astBegin; // start memory management + + // create channel and set encoding + AstFitsChan* chan = astFitsChan(NULL, NULL, ""); + const char* fitswcs = "FITS-WCS"; + const char* encode = currentContext->fits->encoding_; + if (!encode || !*encode) + encode = fitswcs; + astSet (chan, "Card=1, Encoding=%s", encode); + + // ast + AstFrameSet* ast = (AstFrameSet*)astCopy(currentContext->fits->ast_); + + Matrix mx = currentContext->fits->imageToRef * refToWidget * cc; + double ss[] = {mx.matrix(0,0),mx.matrix(1,0), + mx.matrix(0,1),mx.matrix(1,1)}; + double tt[] = {mx.matrix(2,0),mx.matrix(2,1)}; + + AstMatrixMap* mm = astMatrixMap(2, 2, 0, ss, ""); + AstShiftMap* sm = astShiftMap(2, tt, ""); + AstCmpMap* cmp = astCmpMap(mm, sm, 1, ""); + + astRemapFrame(ast, AST__BASE, cmp); + + // write to channel + if (!astWrite(chan, ast)) { + // try again + encode = fitswcs; + astSet (chan, "Card=1, Encoding=%s", encode); + astWrite(chan, ast); } + + // dump cards from channel + astClear(chan, "Card"); + char card[81]; + while (astFindFits(chan, "%f", card, 1)) + dst.cardins(card,NULL); + + astEnd; // now, clean up memory } void FrameBase::saveFitsResampleFits(OutFitsStream& str) -- cgit v0.12