diff options
author | William Joye <wjoye@cfa.harvard.edu> | 2016-10-27 18:59:29 (GMT) |
---|---|---|
committer | William Joye <wjoye@cfa.harvard.edu> | 2016-10-27 18:59:29 (GMT) |
commit | d4d595fa7fb12903db9227d33d48b2b00120dbd1 (patch) | |
tree | 7d18365de0d6d1b29399b6a17c7eb01c2eb3ed49 /tksao/fitsy++/block.C | |
parent | 949f96e29bfe0bd8710d775ce220e597064e2589 (diff) | |
download | blt-d4d595fa7fb12903db9227d33d48b2b00120dbd1.zip blt-d4d595fa7fb12903db9227d33d48b2b00120dbd1.tar.gz blt-d4d595fa7fb12903db9227d33d48b2b00120dbd1.tar.bz2 |
Initial commit
Diffstat (limited to 'tksao/fitsy++/block.C')
-rw-r--r-- | tksao/fitsy++/block.C | 269 |
1 files changed, 269 insertions, 0 deletions
diff --git a/tksao/fitsy++/block.C b/tksao/fitsy++/block.C new file mode 100644 index 0000000..1712b2f --- /dev/null +++ b/tksao/fitsy++/block.C @@ -0,0 +1,269 @@ +// Copyright (C) 1999-2016 +// Smithsonian Astrophysical Observatory, Cambridge, MA, USA +// For conditions of distribution and use, see copyright notice in "copyright" + +#include <iostream> +#include <sstream> +#include <iomanip> +using namespace std; + +#include "block.h" + +FitsBlock::FitsBlock(FitsFile* fits, Vector& block) +{ + FitsHead* srcHead = fits->head(); + + width_ = srcHead->naxis(0)/block[0]; + height_ = srcHead->naxis(1)/block[1]; + if (width_<1) + width_ = 1; + if (height_<1) + height_ = 1; + size_ = (size_t)width_*height_; + + primary_ = fits->primary(); + ext_ = fits->ext(); + inherit_ = fits->inherit(); + + bitpix_ = srcHead->hdu()->bitpix(); + switch (bitpix_) { + case -64: + if (!(data_ = new double[size_])) + return; + dataSize_ = size_*sizeof(double); + break; + default: + bitpix_ = -32; + if (!(data_ = new float[size_])) + return; + dataSize_ = size_*sizeof(float); + break; + } + dataSkip_ = 0; + memset(data_, 0, dataSize_); + + initHeader(fits, block); + + // made it this far, must be valid + byteswap_ = 0; + endian_ = lsb() ? LITTLE : BIG; + valid_ = 1; +} + +FitsBlock::~FitsBlock() +{ + if (data_) + delete [] (char*)data_; +} + +void FitsBlock::initHeader(FitsFile* fits, Vector& block) +{ + head_ = new FitsHead(*(fits->head())); + + // change bitpix + switch (bitpix_) { + case -64: + head_->setInteger("BITPIX", -64, ""); + break; + default: + head_->setInteger("BITPIX", -32, ""); + break; + } + + // change width/height + head_->setInteger("NAXIS1", width_, ""); + head_->setInteger("NAXIS2", height_, ""); + + // IRAF + initLTMV(block); + initCCDSUM(block); + initKeySEC("BIASSEC",block); + initKeySEC("DATASEC",block); + initKeySEC("TRIMSEC",block); + + // WCS + initWCS(block); + + // clear cards + if (head_->find("BZERO")) + head_->carddel("BZERO"); + if (head_->find("BSCALE")) + head_->carddel("BSCALE"); + if (head_->find("DATASUM")) + head_->carddel("DATASUM"); + if (head_->find("CHECKSUM")) + head_->carddel("CHECKSUM"); + if (head_->find("DATAMIN")) + head_->carddel("DATAMIN"); + if (head_->find("DATAMAX")) + head_->carddel("DATAMAX"); + if (head_->find("DATAMEAN")) + head_->carddel("DATAMEAN"); + if (head_->find("GOODMIN")) + head_->carddel("GOODMIN"); + if (head_->find("GOODMAX")) + head_->carddel("GOODMAX"); + if (head_->find("IRAF-MIN")) + head_->carddel("IRAF-MIN"); + if (head_->find("IRAF-MAX")) + head_->carddel("IRAF-MAX"); + + head_->updateHDU(); +} + +void FitsBlock::initCCDSUM(Vector& block) +{ + if (head_->find("CCDSUM")) { + char* val = head_->getString("CCDSUM"); + float xx,yy; + istringstream istr(val); + istr >> xx >> yy; + + xx *= block[0]; + yy *= block[1]; + + ostringstream ostr; + ostr << xx << ' ' << yy << ends; + head_->setString("CCDSUM", ostr.str().c_str(), ""); + } +} + +void FitsBlock::initKeySEC(const char* key, Vector& block) +{ + if (head_->find(key)) { + char* sec = head_->getString(key); + Vector ll,ur; + parseSection(sec,&ll,&ur); + Matrix mm = Translate(-1,-1) * + Scale(1/block[0],1/block[1]) * + Translate(1,1); + + Vector nll = ll*mm; + Vector nur = ur*mm; + ostringstream ostr; + ostr << '[' << int(nll[0]) << ':' << int(nur[0]) + << ',' << int(nll[1]) << ':' << int(nur[1]) << ']' << ends; + head_->setString(key, ostr.str().c_str(), ""); + } +} + +void FitsBlock::initLTMV(Vector& block) +{ + // always do LTMV + double ltv1 = head_->getReal("LTV1", 0); + double ltv2 = head_->getReal("LTV2", 0); + double ltm11 = head_->getReal("LTM1_1", 1); + double ltm12 = head_->getReal("LTM1_2", 0); + double ltm21 = head_->getReal("LTM2_1", 0); + double ltm22 = head_->getReal("LTM2_2", 1); + + Matrix mm(ltm11,ltm12,ltm21,ltm22,ltv1,ltv2); + Matrix im = mm * + Translate(-.5,-.5) * + Scale(1./block[0],1./block[1]) * + Translate(.5,.5); + + if (head_->find("LTV1")) + head_->setReal("LTV1", im.matrix(2,0), 8, ""); + else + head_->appendReal("LTV1", im.matrix(2,0), 8, ""); + + if (head_->find("LTV2")) + head_->setReal("LTV2", im.matrix(2,1), 8, ""); + else + head_->appendReal("LTV2", im.matrix(2,1), 8, ""); + + if (head_->find("LTM1_1")) + head_->setReal("LTM1_1", ltm11/block[0], 8, ""); + else + head_->appendReal("LTM1_1", ltm11/block[0], 8, ""); + + if (head_->find("LTM1_2")) + head_->setReal("LTM1_2", ltm12, 8, ""); + else + head_->appendReal("LTM1_2", ltm12, 8, ""); + + if (head_->find("LTM2_1")) + head_->setReal("LTM2_1", ltm21, 8, ""); + else + head_->appendReal("LTM2_1", ltm21, 8, ""); + + if (head_->find("LTM2_2")) + head_->setReal("LTM2_2", ltm22/block[1], 8, ""); + else + head_->appendReal("LTM2_2", ltm22/block[1], 8, ""); +} + +void FitsBlock::initWCS(Vector& block) +{ + // check for WCS + if (!head_->find("CRPIX1") && !head_->find("CRPIX2")) + return; + + double crpix1 = head_->getReal("CRPIX1", 0); + double crpix2 = head_->getReal("CRPIX2", 0); + Matrix pp = Matrix(1,0,0,1,crpix1,crpix2) * + Translate(-.5,-.5) * + Scale(1/block[0],1/block[1]) * + Translate(.5,.5); + + if (head_->find("CRPIX1")) + head_->setReal("CRPIX1", pp.matrix(2,0), 8, ""); + else + head_->appendReal("CRPIX1", pp.matrix(2,0), 8, ""); + + if (head_->find("CRPIX2")) + head_->setReal("CRPIX2", pp.matrix(2,1), 8, ""); + else + head_->appendReal("CRPIX2", pp.matrix(2,1), 8, ""); + + // CD + if (head_->find("CD1_1") || head_->find("CD1_2") || + head_->find("CD2_1") || head_->find("CD2_2")) { + double cd11 = head_->getReal("CD1_1", 1); + double cd12 = head_->getReal("CD1_2", 0); + double cd21 = head_->getReal("CD2_1", 0); + double cd22 = head_->getReal("CD2_2", 1); + Matrix mm = Matrix(cd11,cd12,cd21,cd22,0,0) * + Scale(block[0],block[1]); + + if (head_->find("CD1_1")) + head_->setReal("CD1_1", mm.matrix(0,0), 8, ""); + else + head_->appendReal("CD1_1", mm.matrix(0,0), 8, ""); + + if (head_->find("CD1_2")) + head_->setReal("CD1_2", mm.matrix(0,1), 8, ""); + else + head_->appendReal("CD1_2", mm.matrix(0,1), 8, ""); + + if (head_->find("CD2_1")) + head_->setReal("CD2_1", mm.matrix(1,0), 8, ""); + else + head_->appendReal("CD2_1", mm.matrix(1,0), 8, ""); + + if (head_->find("CD2_2")) + head_->setReal("CD2_2", mm.matrix(1,1), 8, ""); + else + head_->appendReal("CD2_2", mm.matrix(1,1), 8, ""); + + } + // CDELT + else if (head_->find("CDELT1") || head_->find("CDELT2")) { + double cdelt1 = head_->getReal("CDELT1", 1); + double cdelt2 = head_->getReal("CDELT2", 0); + Matrix mm = Matrix(cdelt1,0,0,cdelt2,0,0) * + Scale(block[0],block[1]); + + if (head_->find("CDELT1")) + head_->setReal("CDELT1", mm.matrix(0,0), 8, ""); + else + head_->appendReal("CDELT1", mm.matrix(0,0), 8, ""); + + if (head_->find("CDELT2")) + head_->setReal("CDELT2", mm.matrix(1,1), 8, ""); + else + head_->appendReal("CDELT2", mm.matrix(1,1), 8, ""); + } +} + |