// Copyright (C) 1999-2018
// Smithsonian Astrophysical Observatory, Cambridge, MA, USA
// For conditions of distribution and use, see copyright notice in "copyright"
#include <pthread.h>
#include "fitsimage.h"
#include "framebase.h"
#include "context.h"
#include "mmap.h"
#include "smmap.h"
#include "mmapincr.h"
#include "alloc.h"
#include "allocgz.h"
#include "channel.h"
#include "share.h"
#include "sshare.h"
#include "socket.h"
#include "socketgz.h"
#include "var.h"
#include "order.h"
#include "iis.h"
#include "hist.h"
#include "compress.h"
#include "analysis.h"
#include "photo.h"
#include "colorscale.h"
#include "wcsast.h"
WCSState::WCSState()
{
wcsSystem_ = Coord::WCS;
wcsSkyFrame_ = Coord::FK5;
};
FitsImage::FitsImage(Context* cx, Tcl_Interp* pp)
{
context_ =cx;
interp_ =pp;
objectKeyword_ =NULL;
fileName =NULL;
rootBaseFileName =NULL;
fullBaseFileName =NULL;
iisFileName =NULL;
fits_ =NULL;
post_ =NULL;
hist_ =NULL;
hpx_ =NULL;
base_ =NULL;
basedata_ =NULL;
manageBlock_ =0;
block_ =NULL;
blockdata_ =NULL;
manageAnalysis_ =0;
analysis_ =NULL;
analysisdata_ =NULL;
image_ =NULL;
data_ =NULL;
nextMosaic_ =NULL;
nextSlice_ =NULL;
keyLTMV =0;
keyATMV =0;
keyDTMV =0;
keyDATASEC =0;
jyBeam_ =1;
imageToData = Translate(-.5, -.5);
dataToImage = Translate( .5, .5);
imageToData3d = Translate3d(-.5, -.5, -.5);
dataToImage3d = Translate3d( .5, .5, .5);
imageToRef3d = imageToData3d;
refToImage3d = dataToImage3d;
manageWCS_ =1;
ast_ =NULL;
encoding_ =NULL;
wcs_ =NULL;
wcsNaxes_ =NULL;
wcsCel_ =NULL;
wcsEqu_ =NULL;
wcsCelLon_ =NULL;
wcsCelLat_ =NULL;
wcsSize_ =NULL;
wcsInv_ =1;
wcsHPX_ =0;
wcsXPH_ =0;
wcsState_ =NULL;
wcsAltHeader_ =NULL;
wfpc2Header_ =NULL;
wcs0Header_ =NULL;
iisMode_ =0;
iiszt_ =0;
for (int ii=0; ii<FTY_MAXAXES; ii++)
address[ii] =1;
}
FitsImage::~FitsImage()
{
if (objectKeyword_)
delete [] objectKeyword_;
if (fileName)
delete [] fileName;
if (rootBaseFileName)
delete [] rootBaseFileName;
if (fullBaseFileName)
delete [] fullBaseFileName;
if (iisFileName)
delete [] iisFileName;
if (fits_)
delete fits_;
if (post_)
delete post_;
if (hist_)
delete hist_;
if (hpx_)
delete hpx_;
if (basedata_)
delete basedata_;
if (manageBlock_) {
if (block_)
delete block_;
if (blockdata_)
delete blockdata_;
}
if (manageAnalysis_) {
if (analysis_)
delete analysis_;
if (analysisdata_)
delete analysisdata_;
}
if (manageWCS_) {
if (ast_)
astAnnul(ast_);
if (encoding_)
delete [] encoding_;
if (wcs_)
delete [] wcs_;
if (wcsNaxes_)
delete [] wcsNaxes_;
if (wcsCel_)
delete [] wcsCel_;
if (wcsEqu_)
delete [] wcsEqu_;
if (wcsCelLon_)
delete [] wcsCelLon_;
if (wcsCelLat_)
delete [] wcsCelLat_;
if (wcsSize_)
delete [] wcsSize_;
if (wcsState_)
delete wcsState_;
}
if (wcsAltHeader_)
delete wcsAltHeader_;
if (wfpc2Header_)
delete wfpc2Header_;
if (wcs0Header_)
delete wcs0Header_;
}
// Fits
FitsImageFitsAlloc::FitsImageFitsAlloc(Context* cx, Tcl_Interp* pp,
const char* ch, const char* fn,
FitsFile::FlushMode flush, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsFitsAlloc(ch, FitsFile::RELAX, flush);
process(fn,id);
}
FitsImageFitsAllocGZ::FitsImageFitsAllocGZ(Context* cx, Tcl_Interp* pp,
const char* ch, const char* fn,
FitsFile::FlushMode flush, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsFitsAllocGZ(ch, FitsFile::RELAX, flush);
process(fn,id);
}
FitsImageFitsChannel::FitsImageFitsChannel(Context* cx, Tcl_Interp* pp,
const char* ch, const char* fn,
FitsFile::FlushMode flush, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsFitsChannel(pp, ch, fn, FitsFile::RELAX, flush);
process(fn,id);
}
FitsImageFitsMMap::FitsImageFitsMMap(Context* cx, Tcl_Interp* pp,
const char* fn, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsFitsMMap(fn, FitsFile::RELAX);
process(fn,id);
}
FitsImageFitsSMMap::FitsImageFitsSMMap(Context* cx, Tcl_Interp* pp,
const char* hdr,
const char* fn, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsFitsSMMap(hdr, fn);
process(fn,id);
}
FitsImageFitsMMapIncr::FitsImageFitsMMapIncr(Context* cx, Tcl_Interp* pp,
const char* fn, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsFitsMMapIncr(fn, FitsFile::RELAX);
process(fn,id);
}
FitsImageFitsShare::FitsImageFitsShare(Context* cx, Tcl_Interp* pp,
Base::ShmType type,
int sid, const char* fn, int id)
: FitsImage(cx, pp)
{
switch (type) {
case Base::SHMID:
fits_ = new FitsFitsShareID(sid, fn, FitsFile::RELAX);
break;
case Base::KEY:
fits_ = new FitsFitsShareKey(sid, fn, FitsFile::RELAX);
break;
}
process(fn,id);
}
FitsImageFitsSShare::FitsImageFitsSShare(Context* cx, Tcl_Interp* pp,
Base::ShmType type,
int hdr, int sid,
const char* fn, int id)
: FitsImage(cx, pp)
{
switch (type) {
case Base::SHMID:
fits_ = new FitsFitsSShareID(hdr, sid, fn);
break;
case Base::KEY:
fits_ = new FitsFitsSShareKey(hdr, sid, fn);
break;
}
process(fn,id);
}
FitsImageFitsSocket::FitsImageFitsSocket(Context* cx, Tcl_Interp* pp,
int s, const char* fn,
FitsFile::FlushMode flush, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsFitsSocket(s, fn, FitsFile::RELAX, flush);
process(fn,id);
}
FitsImageFitsSocketGZ::FitsImageFitsSocketGZ(Context* cx, Tcl_Interp* pp,
int s, const char* fn,
FitsFile::FlushMode flush, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsFitsSocketGZ(s, fn, FitsFile::RELAX, flush);
process(fn,id);
}
FitsImageFitsVar::FitsImageFitsVar(Context* cx, Tcl_Interp* pp,
const char* var, const char* fn, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsFitsVar(pp, var, fn, FitsFile::RELAX);
process(fn,id);
}
FitsImageFitsOrder::FitsImageFitsOrder(Context* cx, Tcl_Interp* pp,
FitsImage* fi,
FitsHead* hdr, char* data, size_t sz,
int id)
: FitsImage(cx, pp)
{
fits_ = new FitsOrder(fi->fitsFile(), hdr, data, sz);
process(NULL,id);
rootBaseFileName = dupstr(fi->rootBaseFileName);
fullBaseFileName = dupstr(fi->fullBaseFileName);
iisFileName = dupstr(fi->fullBaseFileName);
}
// Fits Next
FitsImageFitsNextAlloc::FitsImageFitsNextAlloc(Context* cx, Tcl_Interp* pp,
const char* fn,
FitsFile* prev, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsFitsNextAlloc(prev);
process(fn,id);
}
FitsImageFitsNextAllocGZ::FitsImageFitsNextAllocGZ(Context* cx, Tcl_Interp* pp,
const char* fn,
FitsFile* prev, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsFitsNextAllocGZ(prev);
process(fn,id);
}
FitsImageFitsNextChannel::FitsImageFitsNextChannel(Context* cx, Tcl_Interp* pp,
const char* fn,
FitsFile* prev, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsFitsNextChannel(prev);
process(fn,id);
}
FitsImageFitsNextMMap::FitsImageFitsNextMMap(Context* cx, Tcl_Interp* pp,
const char* fn,
FitsFile* prev, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsFitsNextMMap(prev);
process(fn,id);
}
FitsImageFitsNextSMMap::FitsImageFitsNextSMMap(Context* cx, Tcl_Interp* pp,
const char* fn,
FitsFile* prev, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsFitsNextSMMap(prev);
process(fn,id);
}
FitsImageFitsNextMMapIncr::FitsImageFitsNextMMapIncr(Context* cx, Tcl_Interp* pp,
const char* fn,
FitsFile* prev, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsFitsNextMMapIncr(prev);
process(fn,id);
}
FitsImageFitsNextShare::FitsImageFitsNextShare(Context* cx, Tcl_Interp* pp,
const char* fn,
FitsFile* prev, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsFitsNextShare(prev);
process(fn,id);
}
FitsImageFitsNextSShare::FitsImageFitsNextSShare(Context* cx, Tcl_Interp* pp,
const char* fn,
FitsFile* prev, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsFitsNextSShare(prev);
process(fn,id);
}
FitsImageFitsNextSocket::FitsImageFitsNextSocket(Context* cx, Tcl_Interp* pp,
const char* fn,
FitsFile* prev, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsFitsNextSocket(prev);
process(fn,id);
}
FitsImageFitsNextSocketGZ::FitsImageFitsNextSocketGZ(Context* cx, Tcl_Interp* pp,
const char* fn,
FitsFile* prev, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsFitsNextSocketGZ(prev);
process(fn,id);
}
FitsImageFitsNextVar::FitsImageFitsNextVar(Context* cx, Tcl_Interp* pp,
const char* fn,
FitsFile* prev, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsFitsNextVar(prev);
process(fn,id);
}
FitsImageFitsNextHist::FitsImageFitsNextHist(Context* cx, Tcl_Interp* pp,
FitsImage* fi,
FitsFile* prev, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsHistNext(prev);
process(NULL,id);
fits_->setpFilter(fi->getHistFilter());
fits_->setpBinX(fi->getHistX());
fits_->setpBinY(fi->getHistY());
fits_->setpBinZ(fi->getHistZ());
rootBaseFileName = dupstr(fi->rootBaseFileName);
fullBaseFileName = dupstr(fi->fullBaseFileName);
iisFileName = dupstr(fi->fullBaseFileName);
}
FitsImageFitsNextPost::FitsImageFitsNextPost(Context* cx, Tcl_Interp* pp,
FitsImage* fi,
FitsFile* prev, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsPostNext(prev);
process(NULL,id);
rootBaseFileName = dupstr(fi->rootBaseFileName);
fullBaseFileName = dupstr(fi->fullBaseFileName);
iisFileName = dupstr(fi->fullBaseFileName);
}
FitsImageFitsNextOrder::FitsImageFitsNextOrder(Context* cx, Tcl_Interp* pp,
FitsImage* fi,
FitsFile* prev, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsOrderNext(prev);
process(NULL,id);
rootBaseFileName = dupstr(fi->rootBaseFileName);
fullBaseFileName = dupstr(fi->fullBaseFileName);
iisFileName = dupstr(fi->fullBaseFileName);
}
// Array
FitsImageArrAlloc::FitsImageArrAlloc(Context* cx, Tcl_Interp* pp,
const char* ch, const char* fn,
FitsFile::FlushMode flush, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsArrAlloc(ch, flush);
process(fn,id);
}
FitsImageArrAllocGZ::FitsImageArrAllocGZ(Context* cx, Tcl_Interp* pp,
const char* ch, const char* fn,
FitsFile::FlushMode flush, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsArrAllocGZ(ch, flush);
process(fn,id);
}
FitsImageArrChannel::FitsImageArrChannel(Context* cx, Tcl_Interp* pp,
const char* ch, const char* fn,
FitsFile::FlushMode flush, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsArrChannel(pp, ch, fn, flush);
process(fn,id);
}
FitsImageArrMMap::FitsImageArrMMap(Context* cx, Tcl_Interp* pp,
const char* fn, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsArrMMap(fn);
process(fn,id);
}
FitsImageArrMMapIncr::FitsImageArrMMapIncr(Context* cx, Tcl_Interp* pp,
const char* fn, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsArrMMapIncr(fn);
process(fn,id);
}
FitsImageArrShare::FitsImageArrShare(Context* cx, Tcl_Interp* pp,
Base::ShmType type,
int sid, const char* fn, int id)
: FitsImage(cx, pp)
{
switch (type) {
case Base::SHMID:
fits_ = new FitsArrShareID(sid, fn);
break;
case Base::KEY:
fits_ = new FitsArrShareKey(sid, fn);
break;
}
process(fn,id);
}
FitsImageArrSocket::FitsImageArrSocket(Context* cx, Tcl_Interp* pp,
int s, const char* fn,
FitsFile::FlushMode flush, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsArrSocket(s, fn, flush);
process(fn,id);
}
FitsImageArrSocketGZ::FitsImageArrSocketGZ(Context* cx, Tcl_Interp* pp,
int s, const char* fn,
FitsFile::FlushMode flush, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsArrSocketGZ(s, fn, flush);
process(fn,id);
}
FitsImageArrVar::FitsImageArrVar(Context* cx, Tcl_Interp* pp,
const char* var, const char* fn, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsArrVar(pp, var, fn);
process(fn,id);
}
// ENVI
FitsImageENVISMMap::FitsImageENVISMMap(Context* cx, Tcl_Interp* pp,
const char* hdr,
const char* fn, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsENVISMMap(hdr,fn);
process(fn,id);
}
// NRRD
FitsImageNRRDAlloc::FitsImageNRRDAlloc(Context* cx, Tcl_Interp* pp,
const char* ch, const char* fn,
FitsFile::FlushMode flush, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsNRRDAlloc(ch, flush);
process(fn,id);
}
FitsImageNRRDChannel::FitsImageNRRDChannel(Context* cx, Tcl_Interp* pp,
const char* ch, const char* fn,
FitsFile::FlushMode flush, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsNRRDChannel(pp, ch, fn, flush);
process(fn,id);
}
FitsImageNRRDMMap::FitsImageNRRDMMap(Context* cx, Tcl_Interp* pp,
const char* fn, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsNRRDMMap(fn);
process(fn,id);
}
FitsImageNRRDShare::FitsImageNRRDShare(Context* cx, Tcl_Interp* pp,
Base::ShmType type,
int sid, const char* fn, int id)
: FitsImage(cx, pp)
{
switch (type) {
case Base::SHMID:
fits_ = new FitsNRRDShareID(sid, fn);
break;
case Base::KEY:
fits_ = new FitsNRRDShareKey(sid, fn);
break;
}
process(fn,id);
}
FitsImageNRRDSocket::FitsImageNRRDSocket(Context* cx, Tcl_Interp* pp,
int s, const char* fn,
FitsFile::FlushMode flush, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsNRRDSocket(s, fn, flush);
process(fn,id);
}
FitsImageNRRDVar::FitsImageNRRDVar(Context* cx, Tcl_Interp* pp,
const char* var, const char* fn, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsNRRDVar(pp, var, fn);
process(fn,id);
}
// Photo
FitsImagePhoto::FitsImagePhoto(Context* cx, Tcl_Interp* pp,
const char* ph, const char* fn, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsPhoto(pp, ph);
process(fn,id);
}
FitsImagePhotoCube::FitsImagePhotoCube(Context* cx, Tcl_Interp* pp,
const char* ph, const char* fn, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsPhotoCube(pp, ph);
process(fn,id);
}
FitsImagePhotoCubeNext::FitsImagePhotoCubeNext(Context* cx, Tcl_Interp* pp,
const char* fn,
FitsFile* prev, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsPhotoCubeNext(prev);
process(fn,id);
}
// Mosaic
FitsImageMosaicAlloc::FitsImageMosaicAlloc(Context* cx, Tcl_Interp* pp,
const char* ch,
const char* fn,
FitsFile::FlushMode flush, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsMosaicAlloc(ch, flush);
process(fn,id);
}
FitsImageMosaicAllocGZ::FitsImageMosaicAllocGZ(Context* cx, Tcl_Interp* pp,
const char* ch,
const char* fn,
FitsFile::FlushMode flush,
int id)
: FitsImage(cx, pp)
{
fits_ = new FitsMosaicAllocGZ(ch, flush);
process(fn,id);
}
FitsImageMosaicChannel::FitsImageMosaicChannel(Context* cx, Tcl_Interp* pp,
const char* ch,
const char* fn,
FitsFile::FlushMode flush,
int id)
: FitsImage(cx, pp)
{
fits_ = new FitsMosaicChannel(pp, ch, flush);
process(fn,id);
}
FitsImageMosaicMMap::FitsImageMosaicMMap(Context* cx, Tcl_Interp* pp,
const char* fn, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsMosaicMMap(fn);
process(fn,id);
}
FitsImageMosaicMMapIncr::FitsImageMosaicMMapIncr(Context* cx, Tcl_Interp* pp,
const char* fn,
int id)
: FitsImage(cx, pp)
{
fits_ = new FitsMosaicMMapIncr(fn);
process(fn,id);
}
FitsImageMosaicShare::FitsImageMosaicShare(Context* cx, Tcl_Interp* pp,
Base::ShmType type,
int sid, const char* fn, int id)
: FitsImage(cx, pp)
{
switch (type) {
case Base::SHMID:
fits_ = new FitsMosaicShareID(sid);
break;
case Base::KEY:
fits_ = new FitsMosaicShareKey(sid);
break;
}
process(fn,id);
}
FitsImageMosaicSocket::FitsImageMosaicSocket(Context* cx, Tcl_Interp* pp,
int s, const char* fn,
FitsFile::FlushMode flush, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsMosaicSocket(s, flush);
process(fn,id);
}
FitsImageMosaicSocketGZ::FitsImageMosaicSocketGZ(Context* cx, Tcl_Interp* pp,
int s, const char* fn,
FitsFile::FlushMode flush,
int id)
: FitsImage(cx, pp)
{
fits_ = new FitsMosaicSocketGZ(s, flush);
process(fn,id);
}
FitsImageMosaicVar::FitsImageMosaicVar(Context* cx, Tcl_Interp* pp,
const char* var, const char* fn, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsMosaicVar(pp, var, fn);
process(fn,id);
}
// Mosaic Next
FitsImageMosaicNextAlloc::FitsImageMosaicNextAlloc(Context* cx, Tcl_Interp* pp,
const char* fn,
FitsFile* prev,
FitsFile::FlushMode flush,
int id)
: FitsImage(cx, pp)
{
fits_ = new FitsMosaicNextAlloc(prev, flush);
process(fn,id);
}
FitsImageMosaicNextAllocGZ::FitsImageMosaicNextAllocGZ(Context* cx, Tcl_Interp* pp,
const char* fn,
FitsFile* prev,
FitsFile::FlushMode flush,
int id)
: FitsImage(cx, pp)
{
fits_ = new FitsMosaicNextAllocGZ(prev, flush);
process(fn,id);
}
FitsImageMosaicNextChannel::FitsImageMosaicNextChannel(Context* cx,
Tcl_Interp* pp,
const char* fn,
FitsFile* prev,
FitsFile::FlushMode flush,
int id)
: FitsImage(cx, pp)
{
fits_ = new FitsMosaicNextChannel(prev, flush);
process(fn,id);
}
FitsImageMosaicNextMMap::FitsImageMosaicNextMMap(Context* cx, Tcl_Interp* pp,
const char* fn,
FitsFile* prev, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsMosaicNextMMap(prev);
process(fn,id);
}
FitsImageMosaicNextMMapIncr::FitsImageMosaicNextMMapIncr(Context* cx, Tcl_Interp* pp,
const char* fn,
FitsFile* prev,
int id)
: FitsImage(cx, pp)
{
fits_ = new FitsMosaicNextMMapIncr(prev);
process(fn,id);
}
FitsImageMosaicNextShare::FitsImageMosaicNextShare(Context* cx, Tcl_Interp* pp,
const char* fn,
FitsFile* prev, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsMosaicNextShare(prev);
process(fn,id);
}
FitsImageMosaicNextSocket::FitsImageMosaicNextSocket(Context* cx,
Tcl_Interp* pp,
const char* fn,
FitsFile* prev,
FitsFile::FlushMode flush,
int id)
: FitsImage(cx, pp)
{
fits_ = new FitsMosaicNextSocket(prev, flush);
process(fn,id);
}
FitsImageMosaicNextSocketGZ::FitsImageMosaicNextSocketGZ(Context* cx,
Tcl_Interp* pp,
const char* fn,
FitsFile* prev,
FitsFile::FlushMode flush,
int id)
: FitsImage(cx, pp)
{
fits_ = new FitsMosaicNextSocketGZ(prev, flush);
process(fn,id);
}
FitsImageMosaicNextVar::FitsImageMosaicNextVar(Context* cx, Tcl_Interp* pp,
const char* fn,
FitsFile* prev, int id)
: FitsImage(cx, pp)
{
fits_ = new FitsMosaicNextVar(prev);
process(fn,id);
}
// IIS
FitsImageIIS::FitsImageIIS(Context* cx, Tcl_Interp* pp, int w, int h)
: FitsImage(cx, pp)
{
fits_ = new FitsIIS(w, h);
process("",1);
iisMode_ = 1;
}
void FitsImageIIS::iisErase()
{
((FitsIIS*)fits_)->erase();
}
char* FitsImageIIS::iisGet(int xx, int yy, int dx, int dy)
{
return ((FitsIIS*)fits_)->get(xx, yy, dx, dy);
}
void FitsImageIIS::iisSet(const char* src, int xx, int yy, int dx, int dy)
{
((FitsIIS*)fits_)->set(src, xx, yy, dx, dy);
}
void FitsImageIIS::iisWCS(const Matrix& mm, const Vector& z, int zt)
{
Matrix& mx = (Matrix&)mm;
double sx = mx[0][0];
double sy = mx[1][1];
double tx = mx[2][0];
double ty = mx[2][1];
imageToPhysical =
Translate(0,-height()/2.) *
FlipY() *
Translate(0,height()/2.) *
Translate(-1,0) *
Scale(sx,sy) *
Translate(tx,ty);
physicalToImage = imageToPhysical.invert();
iisz_ = z;
iiszt_ = zt;
}
// FitsImage
void FitsImage::wfpc2WCS(FitsHead* pp, istream& str)
{
FitsHead* hh = parseWCS(str);
// EQUINOX
if (pp->find("EQUINOX")) {
char* equ = pp->getString("EQUINOX");
hh->appendString("EQUINOX", equ, NULL);
}
// DATE-OBS
if (pp->find("DATE-OBS")) {
char* equ = pp->getString("DATE-OBS");
hh->appendString("DATE-OBS", equ, NULL);
}
// Process OBJECT keyword
if (objectKeyword_)
delete [] objectKeyword_;
objectKeyword_ = dupstr(hh->getString("OBJECT"));
// Process WCS keywords
if (wfpc2Header_)
delete wfpc2Header_;
wfpc2Header_ = hh;
initWCS(wfpc2Header_);
}
void FitsImage::appendWCS(istream& str)
{
FitsHead* hh = parseWCS(str);
// process OBJECT keyword
char* obj = dupstr(hh->getString("OBJECT"));
if (obj) {
if (objectKeyword_)
delete [] objectKeyword_;
objectKeyword_ = obj;
}
// Process WCS keywords
FitsHead* hd = image_->head();
// append wcs keywords to the end of the header
int ll = hd->headbytes()+hh->headbytes();
char* cards = new char[ll];
// copy default wcs
memcpy(cards, hd->cards(), hd->headbytes());
// find first END and zero out
for (int i=0; i<hd->headbytes(); i+=80)
if (!strncmp(cards+i,"END",3)) {
memcpy(cards+i, " ",3);
break;
}
// copy appended wcs
memcpy(cards+hd->headbytes(), hh->cards(), hh->headbytes());
delete hh;
if (wcsAltHeader_)
delete wcsAltHeader_;
wcsAltHeader_ = new FitsHead(cards,ll,FitsHead::ALLOC);
initWCS(wcsAltHeader_);
}
char* FitsImage::display(FitsHead* hd)
{
int size = hd->ncard() * (FTY_CARDLEN+1);
char* lbuf = new char[size+1];
char* lptr = lbuf;
char* cptr = hd->cards();
for (int i=0; i<hd->ncard(); i++,cptr+=FTY_CARDLEN) {
memcpy(lptr, cptr, FTY_CARDLEN);
lptr+=FTY_CARDLEN;
*(lptr++) = '\n';
}
lbuf[size] = '\0';
return lbuf;
}
char* FitsImage::displayHeader()
{
Vector blockFactor = context_->blockFactor();
if (blockFactor[0] != 1 && blockFactor[1] != 1)
return display(image_->head());
if (DebugBin || DebugCompress)
return display(image_->head());
else
return display(fits_->head());
}
char* FitsImage::displayPrimary()
{
return display(fits_->primary());
}
char* FitsImage::displayWCS()
{
if (wcsAltHeader_)
return display(wcsAltHeader_);
else if (wfpc2Header_)
return display(wfpc2Header_);
else if (wcs0Header_)
return display(wcs0Header_);
else
return display(image_->head());
}
FitsBound* FitsImage::getDataParams(FrScale::SecMode which)
{
// params in DATA coords 0-n
switch (which) {
case FrScale::IMGSEC:
return &iparams;
case FrScale::DATASEC:
return &dparams;
case FrScale::CROPSEC:
return &cparams;
}
// just to satisfy the compiler
return &iparams;
}
const char* FitsImage::getValue(const Vector& v)
{
if (!isIIS())
return data_->getValue(v);
else {
double val = data_->getValueDouble(v);
ostringstream str;
if (val == 0)
str << ends;
else if (val == IISMIN)
str << '<' << iisz_[0] << ends;
else if (val == IISMAX)
str << '>' << iisz_[1] << ends;
else if (val > IISMAX)
str << ends;
else
// W_LINEAR =1
if (iiszt_ == 1)
str << ((val-IISMIN) * (iisz_[1]-iisz_[0]))/(IISMAX-IISMIN) + iisz_[0]
<< ends;
else
str << val << ends;
memcpy(buf,str.str().c_str(), str.str().length());
return buf;
}
}
void FitsImage::iisSetFileName(const char* fn)
{
if (iisFileName)
delete [] iisFileName;
iisFileName = dupstr(fn);
}
void FitsImage::initWCS(FitsHead* hd)
{
if (manageWCS_) {
if (ast_)
astAnnul(ast_);
ast_ =NULL;
if (encoding_)
delete [] encoding_;
encoding_ =NULL;
if (wcs_)
delete [] wcs_;
wcs_ =NULL;
if (wcsNaxes_)
delete [] wcsNaxes_;
wcsNaxes_ =NULL;
if (wcsCel_)
delete [] wcsCel_;
wcsCel_ =NULL;
if (wcsEqu_)
delete [] wcsEqu_;
wcsEqu_ =NULL;
if (wcsCelLon_)
delete [] wcsCelLon_;
wcsCelLon_ =NULL;
if (wcsCelLat_)
delete [] wcsCelLat_;
wcsCelLat_ =NULL;
if (wcsSize_)
delete [] wcsSize_;
wcsSize_ =NULL;
wcsInv_ = 1;
wcsHPX_ = 0;
wcsXPH_ = 0;
}
// shareWCS?
manageWCS_ =1;
if (context_->shareWCS()) {
FitsImage* ptr = context_->fits;
while (ptr) {
if (ptr == this)
break;
FitsImage* sptr = ptr->nextSlice();
while (sptr) {
if (sptr == this) {
ast_ = ptr->ast_;
encoding_ = ptr->encoding_;
wcs_ = ptr->wcs_;
wcsNaxes_ = ptr->wcsNaxes_;
wcsCel_ = ptr->wcsCel_;
wcsEqu_ = ptr->wcsEqu_;
wcsCelLon_ = ptr->wcsCelLon_;
wcsCelLat_ = ptr->wcsCelLat_;
wcsSize_ = ptr->wcsSize_;
wcsInv_ = ptr->wcsInv_;
wcsHPX_ = ptr->wcsHPX_;
wcsXPH_ = ptr->wcsXPH_;
wcsState_ = ptr->wcsState_;
wcsPhyInit();
manageWCS_ =0;
return;
}
sptr = sptr->nextSlice();
}
ptr = ptr->nextMosaic();
}
}
if (ast_)
astAnnul(ast_);
ast_ =NULL;
if (encoding_)
delete [] encoding_;
encoding_ =NULL;
ast_ = fits2ast(hd);
if (!ast_) {
// reset to process LTMV keywords
keyLTMV =0;
return;
}
// special case
if (astGetI(ast_,"Naxes") == 2 &&
astIsASkyFrame(astGetFrame(ast_,AST__CURRENT)) &&
astGetI(ast_,"LatAxis") == 1) {
int orr[] = {2,1};
astPermAxes(ast_,orr);
}
scanWCS(hd);
// init wcsState
if (wcsState_)
delete wcsState_;
wcsState_ = new WCSState();
astBegin;
wcsSystem(ast_,wcsState_->wcsSystem_);
if (hasWCSEqu(wcsState_->wcsSystem_))
wcsSkyFrame(ast_,wcsState_->wcsSkyFrame_);
astEnd;
// must wait until wcsState_ is realized
wcsSize_ = new double[MULTWCS];
for (int ii=0; ii<MULTWCS; ii++)
wcsSize_[ii] = calcWCSSize((Coord::CoordSystem)(ii+Coord::WCS));
wcsPhyInit();
if (DebugWCS && ast_)
astShow(ast_);
}
void FitsImage::resetWCS()
{
// Process OBJECT keyword
if (objectKeyword_)
delete [] objectKeyword_;
objectKeyword_ = dupstr(image_->getString("OBJECT"));
if (wcsAltHeader_)
delete wcsAltHeader_;
wcsAltHeader_ =NULL;
if (wcs0Header_)
delete wcs0Header_;
wcs0Header_ =NULL;
if (wfpc2Header_)
initWCS(wfpc2Header_);
else
initWCS(image_->head());
// apply block factor
if (ast_) {
Vector block = context_->blockFactor();
if (block[0] != 1 || block[1] != 1) {
astClearStatus; // just to make sure
astBegin; // start memory management
Vector ll(.5,.5);
Vector ur(1.5,1.5);
Vector rr = ur*Translate(-.5,-.5)*Scale(block)*Translate(.5,.5);
AstWinMap* winmap = wcsWinMap(ast_, ll, ur, rr);
if (winmap)
astRemapFrame(ast_, AST__BASE, winmap);
astEnd;
}
}
}
void FitsImage::initWCS0(const Vector& pix)
{
if (!ast_)
return;
FitsHead* hd = new FitsHead(naxis(0), naxis(1), 1, -32);
// CTYPE
hd->appendString("CTYPE1", "RA---TAN", NULL);
hd->appendString("CTYPE2", "DEC--TAN", NULL);
// CRPIX
Vector cc = mapFromRef(pix, Coord::IMAGE, Coord::FK5);
hd->appendReal("CRPIX1", cc[1], 8, NULL);
hd->appendReal("CRPIX2", cc[0], 8, NULL);
// CRVAL
hd->appendReal("CRVAL1", 0, 8, NULL);
hd->appendReal("CRVAL2", 0, 8, NULL);
// CD
double ss = getWCSSize(Coord::WCS);
double ang = getWCSRotation(Coord::WCS,Coord::FK5);
Matrix flip;
switch (getWCSOrientation(Coord::WCS,Coord::FK5)) {
case Coord::NORMAL:
case Coord::YY:
flip = FlipX();
break;
case Coord::XX:
case Coord::XY:
break;
};
Matrix mx = flip*Rotate(ang)*Scale(ss);
hd->appendReal("CD1_1", mx[0][0], 8, NULL);
hd->appendReal("CD1_2", mx[0][1], 8, NULL);
hd->appendReal("CD2_1", mx[1][0], 8, NULL);
hd->appendReal("CD2_2", mx[1][1], 8, NULL);
// EPOCH, EQUINOX
hd->appendReal("EPOCH", 2000, 8, NULL);
hd->appendReal("EQUINOX", 2000, 8, NULL);
// RADESYS
hd->appendString("RADESYS", "FK5", NULL);
if (wcs0Header_)
delete wcs0Header_;
wcs0Header_ = hd;
initWCS(wcs0Header_);
}
void FitsImage::load()
{
if (post_)
base_ = post_;
else if (hpx_)
base_ = hpx_;
else if (hist_)
base_ = hist_;
else
base_ = fits_;
if (basedata_)
delete basedata_;
switch (base_->head()->bitpix()) {
case 8:
basedata_ = new FitsDatam<unsigned char>(base_, interp_);
break;
case 16:
basedata_ = new FitsDatam<short>(base_, interp_);
break;
case -16:
basedata_ = new FitsDatam<unsigned short>(base_, interp_);
break;
case 32:
basedata_ = new FitsDatam<int>(base_, interp_);
break;
case 64:
basedata_ = new FitsDatam<long long>(base_, interp_);
break;
case -32:
basedata_ = new FitsDatam<float>(base_, interp_);
break;
case -64:
basedata_ = new FitsDatam<double>(base_, interp_);
break;
}
// block
block_ = base_;
blockdata_ = basedata_;
// analysis
analysis_ = block_;
analysisdata_ = blockdata_;
// image
image_ = analysis_;
data_ = analysisdata_;
}
void FitsImage::match(const char* xxname1, const char* yyname1,
Coord::CoordSystem sys1, Coord::SkyFrame sky1,
const char* xxname2, const char* yyname2,
Coord::CoordSystem sys2, Coord::SkyFrame sky2,
double rad, Coord::CoordSystem sys,
Coord::DistFormat dist,
const char* rrname)
{
if (!hasWCS(sys1) || !hasWCS(sys2))
return;
// only good for skyframe
astClearStatus; // just to make sure
astBegin; // start memory management
// get lists
Tcl_Obj* listxx1 =
Tcl_GetVar2Ex(interp_, xxname1, NULL, TCL_LEAVE_ERR_MSG);
Tcl_Obj* listyy1 =
Tcl_GetVar2Ex(interp_, yyname1, NULL, TCL_LEAVE_ERR_MSG);
Tcl_Obj* listxx2 =
Tcl_GetVar2Ex(interp_, xxname2, NULL, TCL_LEAVE_ERR_MSG);
Tcl_Obj* listyy2 =
Tcl_GetVar2Ex(interp_, yyname2, NULL, TCL_LEAVE_ERR_MSG);
// get objects
int nxx1;
Tcl_Obj **objxx1;
Tcl_ListObjGetElements(interp_, listxx1, &nxx1, &objxx1);
int nyy1;
Tcl_Obj **objyy1;
Tcl_ListObjGetElements(interp_, listyy1, &nyy1, &objyy1);
int nxx2;
Tcl_Obj **objxx2;
Tcl_ListObjGetElements(interp_, listxx2, &nxx2, &objxx2);
int nyy2;
Tcl_Obj **objyy2;
Tcl_ListObjGetElements(interp_, listyy2, &nyy2, &objyy2);
// sanity check
if (nxx1 != nyy1 || nxx2 != nyy2)
return;
if (!hasWCSCel(sys1) || !hasWCSCel(sys2))
return;
// get doubles
double* ixx1 = new double[nxx1];
for (int ii=0 ; ii<nxx1 ; ii++)
Tcl_GetDoubleFromObj(interp_, objxx1[ii], ixx1+ii);
double* iyy1 = new double[nyy1];
for (int ii=0 ; ii<nyy1 ; ii++)
Tcl_GetDoubleFromObj(interp_, objyy1[ii], iyy1+ii);
double* ixx2 = new double[nxx2];
for (int ii=0 ; ii<nxx2 ; ii++)
Tcl_GetDoubleFromObj(interp_, objxx2[ii], ixx2+ii);
double* iyy2 = new double[nyy2];
for (int ii=0 ; ii<nyy2 ; ii++)
Tcl_GetDoubleFromObj(interp_, objyy2[ii], iyy2+ii);
Vector* in1 = new Vector[nxx1];
for (int ii=0; ii<nxx1; ii++)
in1[ii] = vDegToRad(Vector(ixx1[ii],iyy1[ii]),sys1);
Vector* ptr2 = new Vector[nxx2];
for (int ii=0; ii<nxx2; ii++)
ptr2[ii] = vDegToRad(Vector(ixx2[ii],iyy2[ii]),sys2);
double rr = rad;
switch (dist) {
case Coord::DEGREE:
break;
case Coord::ARCMIN:
rr /= 60.;
break;
case Coord::ARCSEC:
rr /= 60.*60.;
break;
}
rr = zeroTWOPI(degToRad(rr));
Vector* ptr1 =NULL;
if (sky1 != sky2) {
AstFrameSet* wcs1 = (AstFrameSet*)astCopy(ast_);
wcsSystem(wcs1,sys1);
if (hasWCSEqu(sys1))
wcsSkyFrame(wcs1,sky1);
AstFrameSet* wcs2 = (AstFrameSet*)astCopy(ast_);
wcsSystem(wcs2,sys2);
if (hasWCSEqu(sys2))
wcsSkyFrame(wcs2,sky2);
AstFrameSet* cvt = (AstFrameSet*)astConvert(wcs1, wcs2, "SKY");
if (cvt != AST__NULL) {
ptr1 = new Vector[nxx1];
wcsTran(cvt, nxx1, in1, 1, ptr1);
}
}
else
ptr1 = in1;
// now compare
if (ptr1 && ptr2) {
AstFrameSet* wcs = (AstFrameSet*)astCopy(ast_);
wcsSystem(wcs,sys2);
if (hasWCSEqu(sys2))
wcsSkyFrame(wcs,sky2);
Tcl_Obj* objrr = Tcl_NewListObj(0,NULL);
for(int jj=0; jj<nxx2; jj++) {
for (int ii=0; ii<nxx1; ii++) {
double dd = wcsDistance(wcs,ptr1[ii],ptr2[jj]);
if ((dd != AST__BAD) && (dd <= rr)) {
Tcl_Obj* obj[2];
obj[0] = Tcl_NewIntObj(ii+1);
obj[1] = Tcl_NewIntObj(jj+1);
Tcl_Obj* list = Tcl_NewListObj(2,obj);
Tcl_ListObjAppendElement(interp_, objrr, list);
}
}
}
Tcl_SetVar2Ex(interp_, rrname, NULL, objrr, TCL_LEAVE_ERR_MSG);
}
// clean up
astEnd; // now, clean up memory
// clean up
if (ixx1)
delete [] ixx1;
if (iyy1)
delete [] iyy1;
if (ixx2)
delete [] ixx2;
if (iyy2)
delete [] iyy2;
if (!ptr1 && ptr1!=in1)
delete [] ptr1;
if (in1)
delete [] in1;
if (ptr2)
delete [] ptr2;
}
Matrix& FitsImage::matrixToData(Coord::InternalSystem sys)
{
switch (sys) {
case Coord::WINDOW:
return windowToData;
case Coord::CANVAS:
return canvasToData;
case Coord::WIDGET:
return widgetToData;
case Coord::PANNER:
return pannerToData;
case Coord::MAGNIFIER:
return magnifierToData;
case Coord::PS:
return psToData;
case Coord::USER:
return userToData;
case Coord::REF:
return refToData;
}
// just to satisfy the compiler
return refToData;
}
Matrix& FitsImage::matrixFromData(Coord::InternalSystem sys)
{
switch (sys) {
case Coord::WINDOW:
return dataToWindow;
case Coord::CANVAS:
return dataToCanvas;
case Coord::WIDGET:
return dataToWidget;
case Coord::PANNER:
return dataToPanner;
case Coord::MAGNIFIER:
return dataToMagnifier;
case Coord::PS:
return dataToPS;
case Coord::USER:
return dataToUser;
case Coord::REF:
return dataToRef;
}
// just to satisfy the compiler
return dataToRef;
}
Matrix3d& FitsImage::matrixFromData3d(Coord::InternalSystem sys)
{
switch (sys) {
case Coord::WINDOW:
return dataToWindow3d;
case Coord::CANVAS:
return dataToCanvas3d;
case Coord::WIDGET:
return dataToWidget3d;
case Coord::PANNER:
return dataToPanner3d;
case Coord::MAGNIFIER:
return dataToMagnifier3d;
case Coord::PS:
return dataToPS3d;
case Coord::USER:
return dataToUser3d;
case Coord::REF:
return dataToRef3d;
}
// just to satisfy the compiler
return dataToRef3d;
}
Matrix3d& FitsImage::matrixToData3d(Coord::InternalSystem sys)
{
switch (sys) {
case Coord::WINDOW:
return windowToData3d;
case Coord::CANVAS:
return canvasToData3d;
case Coord::WIDGET:
return widgetToData3d;
case Coord::PANNER:
return pannerToData3d;
case Coord::MAGNIFIER:
return magnifierToData3d;
case Coord::PS:
return psToData3d;
case Coord::USER:
return userToData3d;
case Coord::REF:
return refToData3d;
}
// just to satisfy the compiler
return refToData3d;
}
FitsHead* FitsImage::parseWCS(istream& str)
{
FitsHead* hd = image_->head();
FitsHead* rr = new FitsHead(hd->naxis(0), hd->naxis(1),
hd->naxis(2), hd->bitpix());
// this works for both terminated (\n) and non-terminated (fits) headers
while (!str.eof()) {
char buf[256];
str.get(buf,80,'\n');
if (strlen(buf) > 0) {
// check for blank lines
if (*buf == ' ')
break;
string x(buf);
istringstream sstr(x);
char keyword[32];
sstr >> keyword;
if (strchr(buf,'=')) {
char dummy;
sstr >> dummy;
}
if (strchr(buf,'\'')) {
char val[64];
char* ss = strchr(buf,'\'')+1;
char* ee = strrchr(buf,'\'');
int ll = ee-ss;
if (ll<0 || ll>63)
ll =0;
strncpy(val,ss,ll);
val[ll] = '\0';
rr->appendString(keyword, val, "");
}
else {
double val;
sstr >> val;
rr->appendReal(keyword, val, 15, "");
}
if (strlen(buf) <= 80) {
// eat the \n
char b;
str.get(b);
}
}
else
break;
}
return rr;
}
void FitsImage::process(const char* fn, int id)
{
if (!fits_->isValid()) {
reset();
return;
}
if (fits_->isImage()) {
switch (fits_->pEncoding()) {
case FitsFile::RAW:
case FitsFile::BSQ:
break;
case FitsFile::GZIP:
initNRRD();
if (!post_ || !post_->isValid()) {
reset();
return;
}
break;
case FitsFile::BIL:
case FitsFile::BIP:
initENVI();
if (!post_ || !post_->isValid()) {
reset();
return;
}
break;
default:
reset();
return;
}
load();
}
else if (fits_->isBinTable()) {
// Compress
if (fits_->find("ZIMAGE")) {
initCompress();
if (!post_ || !post_->isValid()) {
reset();
return;
}
load();
}
// HEALPIX
else if ((fits_->find("PIXTYPE") && (!strncmp(fits_->getString("PIXTYPE"),"HEALPIX",4))) || fits_->find("NSIDE")) {
initHPX();
if (!hpx_ || !hpx_->isValid()) {
reset();
return;
}
load();
}
else {
// Bintable
initBin();
if (!hist_ || !hist_->isValid()) {
reset();
return;
}
}
}
else if (fits_->isAsciiTable()) {
// HEALPIX
if (fits_->find("NSIDE")) {
initHPX();
if (!hpx_ || !hpx_->isValid()) {
reset();
return;
}
load();
}
}
// set slice address
for (int ii=1; ii<id; ii++) {
for (int jj=2; jj<FTY_MAXAXES; jj++) {
if (address[jj]<naxis(jj)) {
address[jj]++;
break;
}
else
address[jj]=1;
}
}
// load() can call reset()
if (fits_)
setFileName(fn);
}
void FitsImage::processKeywordsPhysical()
{
// Physical to Image (LTM/LTV keywords) (with no wcsname already located)
if (!keyLTMV) {
if (image_->find("LTM1_1") ||
image_->find("LTM1_2") ||
image_->find("LTM2_1") ||
image_->find("LTM2_2") ||
image_->find("LTV1") ||
image_->find("LTV2"))
keyLTMV = 1;
double ltm11 = image_->getReal("LTM1_1", 1);
double ltm12 = image_->getReal("LTM1_2", 0);
double ltm21 = image_->getReal("LTM2_1", 0);
double ltm22 = image_->getReal("LTM2_2", 1);
double ltv1 = image_->getReal("LTV1", 0);
double ltv2 = image_->getReal("LTV2", 0);
physicalToImage = Matrix(ltm11, ltm12, ltm21, ltm22, ltv1, ltv2);
imageToPhysical = physicalToImage.invert();
physicalToImage3d =
Matrix3d(ltm11, ltm12, 0, ltm21, ltm22, 0, 0, 0, 1, ltv1, ltv2, 0);
imageToPhysical3d = physicalToImage3d.invert();
}
// CDD to Detector (DTM/DTV keywords)
keyDTMV =0;
if (image_->find("DTM1_1") ||
image_->find("DTM1_2") ||
image_->find("DTM2_1") ||
image_->find("DTM2_2") ||
image_->find("DTV1") ||
image_->find("DTV2"))
keyDTMV = 1;
double dtm11 = image_->getReal("DTM1_1", 1);
double dtm12 = image_->getReal("DTM1_2", 0);
double dtm21 = image_->getReal("DTM2_1", 0);
double dtm22 = image_->getReal("DTM2_2", 1);
double dtv1 = image_->getReal("DTV1", 0);
double dtv2 = image_->getReal("DTV2", 0);
physicalToDetector = Matrix(dtm11, dtm12, dtm21, dtm22, dtv1, dtv2);
detectorToPhysical = physicalToDetector.invert();
physicalToDetector3d =
Matrix3d(dtm11, dtm12, 0, dtm21, dtm22, 0, 0, 0, 1, dtv1, dtv2, 0);
detectorToPhysical3d = physicalToDetector3d.invert();
// Physical to Amplifier (ATM/ATV keywords)
keyATMV =0;
if (image_->find("ATM1_1") ||
image_->find("ATM1_2") ||
image_->find("ATM2_1") ||
image_->find("ATM2_2") ||
image_->find("ATV1") ||
image_->find("ATV2"))
keyATMV = 1;
double atm11 = image_->getReal("ATM1_1", 1);
double atm12 = image_->getReal("ATM1_2", 0);
double atm21 = image_->getReal("ATM2_1", 0);
double atm22 = image_->getReal("ATM2_2", 1);
double atv1 = image_->getReal("ATV1", 0);
double atv2 = image_->getReal("ATV2", 0);
physicalToAmplifier = Matrix(atm11, atm12, atm21, atm22, atv1, atv2);
amplifierToPhysical = physicalToAmplifier.invert();
physicalToAmplifier3d =
Matrix3d(atm11, atm12, 0, atm21, atm22, 0, 0, 0, 1, atv1, atv2, 0);
amplifierToPhysical3d = physicalToAmplifier3d.invert();
if (DebugMosaic) {
cerr << endl;
cerr << "ATM/V: " << physicalToAmplifier << endl;
cerr << "ATM/V-1: " << amplifierToPhysical << endl;
cerr << "DTM/V: " << physicalToDetector << endl;
cerr << "DTM/V-1: " << detectorToPhysical << endl;
cerr << "LTM/V: " << physicalToImage << endl;
cerr << "LTM/V-1: " << imageToPhysical << endl;
}
/*
// Radio data?
char* bunit = image_->getString("BUNIT");
double cdelt1 = fabs(image_->getReal("CDELT1",0));
double cdelt2 = fabs(image_->getReal("CDELT2",0));
double bmaj = image_->getReal("BMAJ",0);
double bmin = image_->getReal("BMIN",0);
// ok we have a radio map
#define GFACTOR (2.0*sqrt(2.0*log(2.0)))
if (!strncmp(bunit,"JY/BEAM",7) && cdelt1 && cdelt2 && bmaj && bmin) {
// convert from deg to arcsec?
cdelt1 *= 3600;
cdelt2 *= 3600;
bmaj *= 3600;
bmin *= 3600;
jyBeam_ = (2*M_PI*bmaj*bmin)/(GFACTOR*GFACTOR*cdelt1*cdelt2);
}
if (bunit)
delete [] bunit;
if (cunit1)
delete [] cunit1;
if (cunit2)
delete [] cunit2;
*/
}
void FitsImage::processKeywordsParams()
{
// params in DATA coords 0-n
iparams.set(0, 0, width(), height());
{
char* datstr = image_->getString("DATASEC");
// default
Vector v1(1,1);
Vector v2(size());
keyDATASEC =0;
if (datstr && *datstr && parseSection(datstr,&v1,&v2)) {
// additional check
if (v1[0]<1 || v1[1]<1 ||
v1[1]>width() || v2[1]>height() ||
v1[0]>v2[0] || v1[1]>v2[1]) {
// default
v1 = Vector(1,1);
v2 = Vector(size());
keyDATASEC = 0;
} else
keyDATASEC = 1;
}
// dparams is a BBOX in DATA coords 0-n
datasec = BBox(v1,v2);
v1 -= Vector(1,1);
dparams.set(v1[0],v1[1],v2[0],v2[1]);
}
// DEBUG
if (DebugCrop) {
cerr << "iparams " << iparams << endl;
cerr << "dparams " << dparams << endl;
}
}
void FitsImage::processKeywordsFitsSection()
{
Vector ll(dparams.xmin,dparams.ymin);
Vector ur(dparams.xmax,dparams.ymax);
if (fits_->pcoord() && fits_->pxvalid() && fits_->pyvalid()) {
ll[0] = fits_->pxmin();
ur[0] = fits_->pxmax();
ll[1] = fits_->pymin();
ur[1] = fits_->pymax();
ll = ll*physicalToImage*Translate(-1,-1);
ur = ur*physicalToImage;
context_->setSecMode(FrScale::CROPSEC);
}
if (!fits_->pcoord() && fits_->pxvalid()) {
ll[0] = fits_->pxmin()-1;
ur[0] = fits_->pxmax();
context_->setSecMode(FrScale::CROPSEC);
}
if (!fits_->pcoord() && fits_->pyvalid()) {
ll[1] = fits_->pymin()-1;
ur[1] = fits_->pymax();
context_->setSecMode(FrScale::CROPSEC);
}
// params is a BBOX in DATA coords 0-n
setCropParams(ll,ur,0);
// DEBUG
if (DebugCrop)
cerr << "cparams " << cparams << endl;
if (fits_->pzvalid()) {
int zmin = fits_->pzmin()-1;
int zmax = fits_->pzmax();
context_->setSecMode(FrScale::CROPSEC);
context_->setCrop3dParams(zmin,zmax);
}
}
int FitsImage::processKeywordsIRAF(FitsImage* fits)
{
// DETSEC
Coord::Orientation orientation = Coord::NORMAL;
char* detstr = image_->getString("DETSEC");
Vector dv1,dv2;
if (!(detstr && *detstr && parseSection(detstr,&dv1,&dv2)))
return 0;
BBox detsec = BBox(dv1,dv2);
int xx = (dv1[0] < dv2[0]);
int yy = (dv1[1] < dv2[1]);
if (xx && yy)
orientation = Coord::NORMAL;
else if (!xx && yy)
orientation = Coord::XX;
else if (!xx && !yy)
orientation = Coord::XY;
else if (xx && !yy)
orientation = Coord::YY;
// DETSIZE
char* sizestr = image_->getString("DETSIZE");
Vector sv1(1,1);
Vector sv2(10000,10000);
if (sizestr && *sizestr) {
if (!(parseSection(sizestr,&sv1,&sv2)))
return 0;
}
BBox detsize = BBox(sv1,sv2);
// CCDSUM
Vector ccdsum(1,1);
char* ccdstr = image_->getString("CCDSUM");
if (ccdstr && *ccdstr) {
double Ns, Np, Ns1, Np1;
string x(ccdstr);
istringstream str(x);
str >> Ns >> Np >> Ns1 >> Np1;
ccdsum = Vector(1/Ns, 1/Np);
}
// origin
Vector origin = detsec.ll * Scale(ccdsum) * Translate(-datasec.ll);
// matrix
// if the segment is flipped, we can have a discontinuity at
// the edges, due to round off errors, so we 'nudge' it
Matrix flip;
switch (orientation) {
case Coord::NORMAL:
break;
case Coord::XX:
flip = FlipX();
break;
case Coord::YY:
flip = FlipY();
break;
case Coord::XY:
flip = FlipXY();
break;
}
// internal flip
Matrix mflip;
switch (context_->IRAFOrientation(orientation)) {
case Coord::NORMAL:
break;
case Coord::XX:
mflip = FlipX();
break;
case Coord::YY:
mflip = FlipY();
break;
case Coord::XY:
mflip = FlipXY();
break;
}
Vector center = datasec.center() * imageToData;
Vector mcenter = detsize.center() * imageToData * Scale(ccdsum);
wcsToRef_ =
Translate(-center) *
flip *
Translate(center) *
Translate(origin) *
Translate(-mcenter) *
mflip *
Translate(mcenter);
// we do this to shift the origin to the middle of the image to match
// the wcs case. Needed by imageBBox()
// first? reset wcsToRef
if (fits == this) {
irafToRef = wcsToRef_.invert();
wcsToRef_ = Matrix();
}
else
wcsToRef_ *= fits->irafToRef;
if (DebugMosaic) {
cerr << "ProcessKeywordsIRAF" << endl
<< " datasec: " << datasec << endl
<< " ccdsum : " << ccdsum << endl
<< " detsize: " << detsize << endl
<< " detsec : " << detsec << endl
<< " matrix : " << wcsToRef_ << endl;
}
return 1;
}
void FitsImage::replaceWCS(istream& str)
{
FitsHead* hh = parseWCS(str);
// Process OBJECT keyword
if (objectKeyword_)
delete [] objectKeyword_;
objectKeyword_ = dupstr(hh->getString("OBJECT"));
// Process WCS keywords
if (wcsAltHeader_)
delete wcsAltHeader_;
wcsAltHeader_ = hh;
initWCS(wcsAltHeader_);
}
void FitsImage::reset()
{
if (fits_)
delete fits_;
fits_ =NULL;
if (post_)
delete post_;
post_ =NULL;
if (hpx_)
delete hpx_;
hpx_ =NULL;
if (hist_)
delete hist_;
hist_ =NULL;
base_ =NULL;
if (basedata_)
delete basedata_;
basedata_ =NULL;
if (manageBlock_) {
if (block_)
delete block_;
if (blockdata_)
delete blockdata_;
}
manageBlock_ =0;
block_ =NULL;
blockdata_ =NULL;
if (manageAnalysis_) {
if (analysis_)
delete analysis_;
if (analysisdata_)
delete analysisdata_;
}
manageAnalysis_ =0;
analysis_ =NULL;
analysisdata_ =NULL;
image_ =NULL;
data_ =NULL;
}
char* FitsImage::root(const char* fn)
{
if (fn) {
const char* ptr = fn; // init the ptr
while(*ptr++); // walk it forward to end of string
ptr--; // backup one
while(*ptr != '/' && ptr != fn) // walk it backward til last / or beginning
ptr--;
if (*ptr == '/') // step it over the last '/'
ptr++;
return dupstr(ptr); // got it!
}
else
return NULL;
}
void FitsImage::setCropParams(int datasec)
{
// params in DATA coords 0-n
if (!datasec)
cparams = iparams;
else
cparams = dparams;
}
void FitsImage::setCropParams(const Vector& ss, const Vector& tt, int datasec)
{
// Coord are in DATA
Vector ll = ss;
Vector ur = tt;
int xmin = ll[0];
int xmax = ur[0];
int ymin = ll[1];
int ymax = ur[1];
if (xmin>xmax) {
xmin = ur[0];
xmax = ll[0];
}
if (ymin>ymax) {
ymin = ur[1];
ymax = ll[1];
}
setCropParams(xmin,ymin,xmax,ymax,datasec);
}
void FitsImage::setCropParams(int x0, int y0, int x1, int y1, int datasec)
{
// params in DATA coords 0-n
FitsBound* params;
if (!datasec)
params = &iparams;
else
params = &dparams;
if (x0<params->xmin)
x0=params->xmin;
if (x0>params->xmax)
x0=params->xmax;
if (x1<params->xmin)
x1=params->xmin;
if (x1>params->xmax)
x1=params->xmax;
if (y0<params->ymin)
y0=params->ymin;
if (y0>params->ymax)
y0=params->ymax;
if (y1<params->ymin)
y1=params->ymin;
if (y1>params->ymax)
y1=params->ymax;
cparams.set(x0,y0,x1,y1);
}
void FitsImage::setFileName(const char* fn)
{
if (fileName)
delete [] fileName;
fileName = NULL;
if (rootBaseFileName)
delete [] rootBaseFileName;
rootBaseFileName = NULL;
if (fullBaseFileName)
delete [] fullBaseFileName;
fullBaseFileName = NULL;
if (iisFileName)
delete [] iisFileName;
iisFileName = NULL;
// no filename to set
if (!fn)
return;
// strip any '[]'
char* ffn = strip(fn);
FitsFile* ptr = post_ ? post_ : fits_;
if (!ptr)
return;
const char* ext = ptr->extname();
if (ext) {
{
ostringstream str;
str << ffn << '[' << ext << ']' << ends;
fullBaseFileName = dupstr(str.str().c_str());
}
{
char* m = root(ffn);
ostringstream str;
str << m << '[' << ext << ']' << ends;
rootBaseFileName = dupstr(str.str().c_str());
delete [] m;
}
}
else if (ptr->ext()) {
{
ostringstream str;
str << ffn << '[' << ptr->ext() << ']' << ends;
fullBaseFileName = dupstr(str.str().c_str());
}
{
char* m = root(ffn);
ostringstream str;
str << m << '[' << ptr->ext() << ']' << ends;
rootBaseFileName = dupstr(str.str().c_str());
delete [] m;
}
}
else {
fullBaseFileName = dupstr(ffn);
rootBaseFileName = root(ffn);
}
// by default, iisFileName is fullBaseFileName
if (fullBaseFileName)
iisFileName = dupstr(fullBaseFileName);
delete [] ffn;
}
void FitsImage::setObjectKeyword(const char* obj)
{
if (objectKeyword_)
delete [] objectKeyword_;
objectKeyword_ = dupstr(obj);
}
char* FitsImage::strip(const char* fn)
{
if (fn) {
char* r = dupstr(fn); // dup the string
char* ptr = r; // init the ptr
while(*ptr != '[' && *ptr) // walk it forward til '[' or end
ptr++;
*ptr = '\0'; // zero out rest
return r; // got it!
}
else
return NULL;
}
int FitsImage::nhdu()
{
int dd =1;
for (int ii=2; ii<FTY_MAXAXES; ii++)
if (naxis(ii))
dd *= naxis(ii);
return dd;
}
void FitsImage::updateClip(FrScale* fr)
{
data_->updateClip(fr,getDataParams(fr->secMode()));
}
void* clipproc(void* tt)
{
t_clip_arg* targ = (t_clip_arg*)tt;
FitsData* data = targ->data;
FrScale* fr = targ->fr;
FitsBound* bb = targ->bb;
data->updateClip(fr,bb);
return NULL;
}
void FitsImage::updateClip(FrScale* fr, pthread_t* thread, t_clip_arg* targ)
{
targ->data = data_;
targ->fr = fr;
targ->bb = getDataParams(fr->secMode());
int result = pthread_create(thread, NULL, clipproc, targ);
if (result)
internalError("Unable to Create Thread");
}
const char* FitsImage::getFileName(Base::FileNameType type)
{
switch (type) {
case Base::ROOTBASE:
return rootBaseFileName;
case Base::FULLBASE:
return fullBaseFileName;
case Base::ROOT:
case Base::FULL:
// clear the buffer
if (fileName)
delete [] fileName;
fileName =NULL;
// generate filename
// if FITS bin table cube, be sure to check the first slice
if (context_->fits->isHist())
return updateFileNameBin(type);
else
return updateFileNameImage(type);
}
// just to satisfy the compiler
return rootBaseFileName;
}
const char* FitsImage::updateFileNameImage(Base::FileNameType type)
{
// 2d/3d section
char* sec =NULL;
switch (context_->secMode()) {
case FrScale::IMGSEC:
case FrScale::DATASEC:
{
Vector blockFactor = context_->blockFactor();
if (blockFactor[0] != 1) {
ostringstream str;
str << "*," << blockFactor[0] << ends;
sec = dupstr(str.str().c_str());
}
}
break;
case FrScale::CROPSEC:
{
FitsBound* params =getDataParams(FrScale::CROPSEC);
// params is a BBOX in DATA coords 0-n
// xlate to 1-n
Vector ll= Vector(params->xmin,params->ymin) * Translate(1,1);
Vector ur(params->xmax,params->ymax);
Vector blockFactor = context_->blockFactor();
if (blockFactor[0] != 1) {
ostringstream str;
str << ll[0] << ':' << ur[0] << ','
<< ll[1] << ':' << ur[1] << ','
<< blockFactor[0] << ends;
sec = dupstr(str.str().c_str());
}
else {
ostringstream str;
str << ll[0] << ':' << ur[0] << ','
<< ll[1] << ':' << ur[1] << ends;
sec = dupstr(str.str().c_str());
}
}
break;
}
// address
char* add =NULL;
{
int doAdd =0;
ostringstream str;
int jj;
for (jj=FTY_MAXAXES-1; jj>=2; jj--) {
if (address[jj]!=1)
break;
}
jj++;
for (int ii=2; ii<jj; ii++) {
doAdd =1;
if (ii==2)
str << "plane=";
else
str << ':';
str << address[ii];
}
if (doAdd) {
str << ends;
add = dupstr(str.str().c_str());
}
}
switch (type) {
case Base::ROOTBASE:
case Base::FULLBASE:
// better not get here
return NULL;
case Base::ROOT:
if (rootBaseFileName) {
ostringstream str;
if (add && sec)
str << rootBaseFileName << '[' << add << ']' << '[' << sec << ']';
else if (add && !sec)
str << rootBaseFileName << '[' << add << ']';
else if (!add && sec)
str << rootBaseFileName << '[' << sec << ']';
else
str << rootBaseFileName;
str << ends;
fileName = dupstr(str.str().c_str());
}
break;
case Base::FULL:
if (fullBaseFileName) {
ostringstream str;
if (add && sec)
str << fullBaseFileName << '[' << add << ']' << '[' << sec << ']';
else if (add && !sec)
str << fullBaseFileName << '[' << add << ']';
else if (!add && sec)
str << fullBaseFileName << '[' << sec << ']';
else
str << fullBaseFileName;
str << ends;
fileName = dupstr(str.str().c_str());
}
break;
}
if (sec)
delete [] sec;
if (add)
delete [] add;
return fileName;
}
const char* FitsImage::updateFileNameBin(Base::FileNameType type)
{
char* filter = (char*)fits_->pFilter();
int doFilter = (filter && *filter);
// x,y filter
char* sec = NULL;
{
switch (context_->secMode()) {
case FrScale::IMGSEC:
case FrScale::DATASEC:
break;
case FrScale::CROPSEC:
{
ostringstream str;
FitsBound* params =getDataParams(FrScale::CROPSEC);
// dataToPhysical not set yet
Vector ll = Vector(params->xmin,params->ymin) *
dataToImage * imageToPhysical;
Vector ur = Vector(params->xmax,params->ymax) *
dataToImage * imageToPhysical;
str << fits_->pBinX() << ">=" << ll[0] << ','
<< fits_->pBinX() << "<=" << ur[0] << ','
<< fits_->pBinY() << ">=" << ll[1] << ','
<< fits_->pBinY() << "<=" << ur[1] << ends;
sec = dupstr(str.str().c_str());
}
break;
}
}
// z filter
char* slice =NULL;
{
char* zcol = (char*)fits_->pBinZ();
int bd = context_->binDepth();
if (bd>1 && zcol) {
// only the first slice will have this
FitsImage* first = context_->fits;
if (first) {
Vector zlim = first->fits_->getColMinMax(zcol);
double zlen = zlim[1]-zlim[0];
double zdelta = zlen/bd;
double zptr = zlim[0] + (address[2]-1)*zdelta;
ostringstream str;
str << zcol << ">=" << zptr << '&' << zcol << '<' << zptr+zdelta <<ends;
slice = dupstr(str.str().c_str());
}
}
}
switch (type) {
case Base::ROOTBASE:
case Base::FULLBASE:
// better not get here
return NULL;
case Base::ROOT:
if (rootBaseFileName) {
ostringstream str;
str << rootBaseFileName;
if (doFilter && sec && slice)
str << '[' << filter << ',' << sec << ',' << slice << ']';
else if (doFilter && sec && !slice)
str << '[' << filter << ',' << sec << ']';
else if (doFilter && !sec && slice)
str << '[' << filter << ',' << slice << ']';
else if (doFilter && !sec && !slice)
str << '[' << filter << ']';
else if (!doFilter && sec && slice)
str << '[' << sec << ',' << slice << ']';
else if (!doFilter && sec && !slice)
str << '[' << sec << ']';
else if (!doFilter && !sec && slice)
str << '[' << slice << ']';
str << ends;
fileName = dupstr(str.str().c_str());
}
case Base::FULL:
if (fullBaseFileName) {
ostringstream str;
str << fullBaseFileName;
if (doFilter && sec && slice)
str << '[' << filter << ',' << sec << ',' << slice << ']';
else if (doFilter && sec && !slice)
str << '[' << filter << ',' << sec << ']';
else if (doFilter && !sec && slice)
str << '[' << filter << ',' << slice << ']';
else if (doFilter && !sec && !slice)
str << '[' << filter << ']';
else if (!doFilter && sec && slice)
str << '[' << sec << ',' << slice << ']';
else if (!doFilter && sec && !slice)
str << '[' << sec << ']';
else if (!doFilter && !sec && slice)
str << '[' << slice << ']';
str << ends;
fileName = dupstr(str.str().c_str());
}
break;
}
return fileName;
}
void FitsImage::updateMatrices(Matrix& rgbToRef,
Matrix& refToUser,
Matrix& userToWidget,
Matrix& widgetToCanvas,
Matrix& canvasToWindow)
{
dataToRef = wcsToRef_ * rgbToRef;
refToData = dataToRef.invert();
dataToUser = dataToRef * refToUser;
userToData = dataToUser.invert();
dataToWidget = dataToUser * userToWidget;
widgetToData = dataToWidget.invert();
dataToCanvas = dataToWidget * widgetToCanvas;
canvasToData = dataToCanvas.invert();
dataToWindow = dataToCanvas * canvasToWindow;
windowToData = dataToWindow.invert();
refToCanvas = refToUser * userToWidget * widgetToCanvas;
canvasToRef = refToCanvas.invert();
imageToRef = imageToData * dataToRef;
refToImage = imageToRef.invert();
imageToWidget = imageToRef * refToUser * userToWidget;
widgetToImage = imageToWidget.invert();
physicalToRef = physicalToImage * imageToData * dataToRef;
refToPhysical = physicalToRef.invert();
amplifierToRef = amplifierToPhysical * physicalToRef;
refToAmplifier = amplifierToRef.invert();
detectorToRef = detectorToPhysical * physicalToRef;
refToDetector = detectorToRef.invert();
physicalToRef3d = physicalToImage3d * imageToData3d * dataToRef3d;
refToPhysical3d = physicalToRef3d.invert();
amplifierToRef3d = amplifierToPhysical3d * physicalToRef3d;
refToAmplifier3d = amplifierToRef3d.invert();
detectorToRef3d = detectorToPhysical3d * physicalToRef3d;
refToDetector3d = detectorToRef3d.invert();
}
void FitsImage::updateMatrices(Matrix3d& refToUser3d,
Matrix3d& userToWidget3d,
Matrix3d& widgetToCanvas3d,
Matrix3d& canvasToWindow3d)
{
dataToUser3d = dataToRef3d * refToUser3d;
userToData3d = dataToUser3d.invert();
dataToWidget3d = dataToUser3d * userToWidget3d;
widgetToData3d = dataToWidget3d.invert();
dataToCanvas3d = dataToWidget3d * widgetToCanvas3d;
canvasToData3d = dataToCanvas3d.invert();
dataToWindow3d = dataToCanvas3d * canvasToWindow3d;
windowToData3d = dataToWindow3d.invert();
}
void FitsImage::updatePannerMatrices(Matrix& refToPanner)
{
dataToPanner = dataToRef * refToPanner;
pannerToData = dataToPanner.invert();
}
void FitsImage::updatePannerMatrices(Matrix3d& refToPanner3d)
{
dataToPanner3d = dataToRef3d * refToPanner3d;
pannerToData3d = dataToPanner3d.invert();
}
void FitsImage::updateMagnifierMatrices(Matrix& refToMagnifier)
{
dataToMagnifier = dataToRef * refToMagnifier;
magnifierToData = dataToMagnifier.invert();
}
void FitsImage::updateMagnifierMatrices(Matrix3d& refToMagnifier3d)
{
dataToMagnifier3d = dataToRef3d * refToMagnifier3d;
magnifierToData3d = dataToMagnifier3d.invert();
}
void FitsImage::updatePS(Matrix ps)
{
dataToPS = dataToRef * ps;
psToData = dataToPS.invert();
}
void FitsImage::updatePS(Matrix3d ps)
{
dataToPS3d = dataToRef3d * ps;
psToData3d = dataToPS3d.invert();
}
// WCS
double FitsImage::getWCSSize(Coord::CoordSystem sys)
{
if (!wcsSize_ || sys<Coord::WCS)
return 0;
else
return wcsSize_[sys-Coord::WCS];
}
double FitsImage::calcWCSSize(Coord::CoordSystem sys)
{
if (!hasWCS(sys))
return 0;
astClearStatus; // just to make sure
setWCSSystem(sys);
// only check lon axis as it will not change near poles
Vector in[3];
Vector out[3];
in[0] = center();
in[1] = center()+Vector(0,1);
wcsTran(ast_, 2, in, 1, out);
double dd = wcsDistance(ast_, out[0], out[1]);
return hasWCSCel(sys) ? radToDeg(dd) : dd;
}
Coord::Orientation FitsImage::getWCSOrientation(Coord::CoordSystem sys,
Coord::SkyFrame sky)
{
if (!hasWCS(sys))
return Coord::NORMAL;
// special case, makes no sense
if (wcsXPH_)
return Coord::NORMAL;
astClearStatus; // just to make sure
astBegin; // start memory management
setWCSSysSkyFrame(sys,sky);
Vector in[3];
Vector out[3];
in[0] = center();
in[1] = center()+Vector(0,1);
in[2] = center()+Vector(1,0);
wcsTran(ast_, 3, in, 1, out);
double ang = wcsAngle(ast_,out[1],out[0],out[2]);
astEnd;
if (!(isnan(ang)||isinf(ang)||(ang == -DBL_MAX)||(ang == DBL_MAX))) {
if ((hasWCSCel(sys) && ang>0) || (!hasWCSCel(sys) && ang<0))
return Coord::XX;
else
return Coord::NORMAL;
}
else
return Coord::NORMAL;
}
double FitsImage::getWCSRotation(Coord::CoordSystem sys, Coord::SkyFrame sky)
{
if (!hasWCS(sys))
return 0;
// special case, makes no sense
if (wcsXPH_)
return 0;
astClearStatus; // just to make sure
astBegin; // start memory management
setWCSSysSkyFrame(sys,sky);
if (!wcsHPX_) {
Vector in[3];
Vector out[3];
in[0] = center();
in[1] = center()+Vector(0,1);
in[2] = center()+Vector(1,0);
wcsTran(ast_, 3, in, 1, out);
double ang = wcsAxAngle(ast_,out[0], out[1]);
double ang3 = wcsAngle(ast_,out[1],out[0],out[2]);
astEnd;
if ((!(isnan(ang )||isinf(ang) ||(ang ==-DBL_MAX)||(ang ==DBL_MAX))) &&
(!(isnan(ang3)||isinf(ang3)||(ang3==-DBL_MAX)||(ang3==DBL_MAX)))) {
if ((hasWCSCel(sys) && ang3>0) || (!hasWCSCel(sys) && ang3<0))
return -ang;
else
return ang;
}
else
return 0;
}
else { // special case for HPX
Vector cc = center();
Vector wcc = wcsTran(ast_, cc, 1);
Vector wnorth = wcc + Vector(0,.001);
Vector north = wcsTran(ast_, wnorth,0);
int current = astGetI(ast_,"Current");
int base = astGetI(ast_,"Base");
astSetI(ast_,"Current",base);
double ang = wcsAxAngle(ast_,cc,north);
astSetI(ast_,"Current",current);
astEnd;
if (!(isnan(ang)||isinf(ang)||(ang == -DBL_MAX)||(ang == DBL_MAX)))
return ang;
else
return 0;
}
}
double FitsImage::getWCSDist(const Vector& vv1, const Vector& vv2,
Coord::CoordSystem sys)
{
if (!hasWCS(sys))
return 0;
astClearStatus; // just to make sure
astBegin; // start memory management
setWCSSystem(sys);
astEnd;
return hasWCSCel(sys) ?
radToDeg(wcsDistance(ast_,vDegToRad(vv1,sys),vDegToRad(vv2,sys))) :
wcsDistance(ast_,vv1,vv2);
}
const char* FitsImage::getWCSDomain(Coord::CoordSystem sys)
{
if (!hasWCS(sys))
return NULL;
astClearStatus; // just to make sure
astBegin; // start memory management
setWCSSystem(sys);
astEnd;
const char* domain = astGetC(ast_, "Domain");
return domain;
}
const char* FitsImage::getWCSSystem(Coord::CoordSystem sys)
{
if (!hasWCSCel(sys))
return NULL;
astClearStatus; // just to make sure
astBegin; // start memory management
AstFrameSet* fs =
(AstFrameSet*)astFindFrame(ast_, astSkyFrame(" MaxAxes=4")," ");
const char* str =NULL;
if (fs)
str = astGetC(fs, "System");
astEnd;
return str;
}
const char* FitsImage::getWCSAxisSymbol(Coord::CoordSystem sys, int axis)
{
if (!hasWCS(sys))
return NULL;
int id = sys-Coord::WCS;
if (wcsNaxes_[id] <= axis)
return NULL;
ostringstream str;
str << "Symbol(" << axis+1 << ")" << ends;
return astGetC(ast_, str.str().c_str());
}
Vector FitsImage::pix2wcs(const Vector& in, Coord::CoordSystem sys,
Coord::SkyFrame sky)
{
if (!hasWCS(sys))
return Vector();
astClearStatus; // just to make sure
astBegin; // start memory management
setWCSSysSkyFrame(sys,sky);
Vector out = wcsTran(ast_, in, 1);
if (!astOK || !checkWCS(out))
return Vector();
astNorm(ast_, out.v);
astEnd;
return vRadToDeg(out,sys);
}
VectorStr FitsImage::pix2wcs(const Vector& in, Coord::CoordSystem sys,
Coord::SkyFrame sky, Coord::SkyFormat format)
{
if (!hasWCS(sys))
return VectorStr();
astClearStatus; // just to make sure
astBegin; // start memory management
setWCSSysSkyFrame(sys,sky);
Vector out = wcsTran(ast_, in, 1);
if (!astOK || !checkWCS(out))
return VectorStr();
setWCSFormat(sys,sky,format);
astNorm(ast_, out.v);
astEnd;
int naxes = astGetI(ast_,"Naxes");
switch (naxes) {
case 1:
{
ostringstream str;
str << setprecision(context_->parent_->precLinear_) << in[1] << ends;
return VectorStr(astFormat(ast_,1,out[0]), str.str().c_str());
}
case 2:
case 3:
case 4:
return VectorStr(astFormat(ast_,1,out[0]), astFormat(ast_,2,out[1]));
}
return VectorStr();
}
Vector3d FitsImage::pix2wcs(const Vector3d& in, Coord::CoordSystem sys,
Coord::SkyFrame sky)
{
if (!(hasWCS(sys) && hasWCS3D(sys)))
return Vector();
astClearStatus; // just to make sure
astBegin; // start memory management
setWCSSysSkyFrame(sys,sky);
Vector3d out = wcsTran(ast_, in, 1);
if (!astOK || !checkWCS(out))
return Vector3d();
astNorm(ast_, out.v);
astEnd;
return vRadToDeg(out,sys);
}
VectorStr3d FitsImage::pix2wcs(const Vector3d& in, Coord::CoordSystem sys,
Coord::SkyFrame sky, Coord::SkyFormat format)
{
if (!hasWCS(sys))
return VectorStr3d();
astClearStatus; // just to make sure
astBegin; // start memory management
setWCSSysSkyFrame(sys,sky);
Vector3d out = wcsTran(ast_, in, 1);
if (!astOK || !checkWCS(out))
return VectorStr3d();
setWCSFormat(sys,sky,format);
astNorm(ast_, out.v);
astEnd;
int naxes = astGetI(ast_,"Naxes");
switch (naxes) {
case 1:
{
ostringstream str1;
str1 << setprecision(context_->parent_->precLinear_) << in[1] << ends;
ostringstream str2;
str2 << setprecision(context_->parent_->precLinear_) << in[2] << ends;
return VectorStr3d(astFormat(ast_,1,out[0]),
str1.str().c_str(),
str2.str().c_str());
}
case 2:
{
ostringstream str;
str << setprecision(context_->parent_->precLinear_) << in[1] << ends;
return VectorStr3d(astFormat(ast_,1,out[0]),
astFormat(ast_,2,out[1]),
str.str().c_str());
}
case 3:
case 4:
return VectorStr3d(astFormat(ast_,1,out[0]),
astFormat(ast_,2,out[1]),
astFormat(ast_,3,out[2]));
}
return VectorStr3d();
}
Vector FitsImage::wcs2pix(const Vector& vv, Coord::CoordSystem sys,
Coord::SkyFrame sky)
{
if (!hasWCS(sys) || !wcsInv_) {
maperr =1;
return Vector();
}
astClearStatus; // just to make sure
astBegin; // start memory management
setWCSSysSkyFrame(sys,sky);
Vector in = vDegToRad(vv,sys);
Vector out = wcsTran(ast_, in, 0);
astEnd;
if (!astOK || !checkWCS(out)) {
maperr =1;
return Vector();
}
maperr =0;
return out;
}
Vector3d FitsImage::wcs2pix(const Vector3d& vv, Coord::CoordSystem sys,
Coord::SkyFrame sky)
{
if (!hasWCS(sys) || !hasWCS3D(sys) || !wcsInv_)
return Vector3d();
astClearStatus; // just to make sure
astBegin; // start memory management
setWCSSysSkyFrame(sys,sky);
Vector3d in = vDegToRad(vv,sys);
Vector3d out = wcsTran(ast_, in, 0);
astEnd;
if (!astOK || !checkWCS(out))
return Vector3d();
return out;
}
int FitsImage::hasWCS(Coord::CoordSystem sys)
{
if (!wcs_ || sys<Coord::WCS)
return 0;
else
return wcs_[sys-Coord::WCS];
}
int FitsImage::hasWCSCel(Coord::CoordSystem sys)
{
if (!wcsCel_ || sys<Coord::WCS)
return 0;
else
return wcsCel_[sys-Coord::WCS];
}
int FitsImage::hasWCSEqu(Coord::CoordSystem sys)
{
if (!wcsEqu_ || sys<Coord::WCS)
return 0;
else
return wcsEqu_[sys-Coord::WCS];
}
int FitsImage::hasWCSLinear(Coord::CoordSystem sys)
{
if (!wcs_ || !wcsCel_ || sys<Coord::WCS)
return 0;
else
return wcs_[sys-Coord::WCS] && !wcsCel_[sys-Coord::WCS];
}
int FitsImage::hasWCS3D(Coord::CoordSystem sys)
{
if (!wcsNaxes_ || sys<Coord::WCS)
return 0;
else
return (wcsNaxes_[sys-Coord::WCS]>2) ? 1 : 0;
}
void FitsImage::scanWCS(FitsHead* hd)
{
// init wcs_ array
wcs_ = new int[MULTWCS];
for (int ii=0; ii<MULTWCS; ii++)
wcs_[ii] =0;
wcsNaxes_ = new int[MULTWCS];
for (int ii=0; ii<MULTWCS; ii++)
wcsNaxes_[ii] =0;
// init wcsCel_ array
wcsCel_ = new int[MULTWCS];
for (int ii=0; ii<MULTWCS; ii++)
wcsCel_[ii] =0;
wcsEqu_ = new int[MULTWCS];
for (int ii=0; ii<MULTWCS; ii++)
wcsEqu_[ii] =0;
wcsCelLon_ = new int[MULTWCS];
for (int ii=0; ii<MULTWCS; ii++)
wcsCelLon_[ii] =0;
wcsCelLat_ = new int[MULTWCS];
for (int ii=0; ii<MULTWCS; ii++)
wcsCelLat_[ii] =0;
if (!ast_)
return;
// easy one, HPX?
char key[] = "CTYPE1 ";
if (image_) {
const char* str = image_->getKeyword(key);
if (str) {
if (!strncmp(str+5,"HPX",3))
wcsHPX_ =1;
if (!strncmp(str+5,"XPH",3))
wcsXPH_ =1;
delete [] str;
}
}
astClearStatus;
astBegin;
int nframes = astGetI(ast_, "Nframe");
// wcs_[]
// since we have ast_
wcs_[0] =1;
wcsNaxes_[0] = astGetI(ast_,"Naxes");
// fill out wcs_ array
for (int kk=0; kk<nframes; kk++) {
AstFrameSet* ff = (AstFrameSet*)astGetFrame(ast_,kk+1);
const char* id = astGetC(ff, "Ident");
if (id && *id) {
int jj = (*id == ' ') ? 0 : *id-'@';
wcs_[jj] = 1;
wcsNaxes_[jj] = astGetI(ff,"Naxes");
}
}
// wcsCel_[]
for (int kk=0; kk<nframes; kk++) {
AstFrame* ff = (AstFrame*)astGetFrame(ast_,kk+1);
int ii=0; // current WCS
const char* id = astGetC(ff, "Ident");
if (id && *id)
ii = (*id == ' ') ? 0 : *id-'@';
AstFrameSet* fs =
(AstFrameSet*)astFindFrame(ff, astSkyFrame(" MaxAxes=4")," ");
if (fs) {
wcsCel_[ii] =1;
const char* sys = astGetC(ff, "System");
if (sys && *sys) {
// cerr << "System: " << sys << endl;
if (strncmp(sys,"Unknown",7))
wcsEqu_[ii] =1;
for (int jj=0; jj<wcsNaxes_[ii]; jj++) {
ostringstream str;
str << "Symbol(" << jj+1 << ")" << ends;
const char* ss = astGetC(ff, str.str().c_str());
size_t lss = strlen(ss);
// cerr << "Symbol: " << ss << endl;
if (ss) {
if (!strcmp(ss,"RA") ||
!strcmp(ss,"l") ||
!strcmp(ss,"Lambda") ||
(lss>1 && !strcmp(ss+1,"LON")) ||
(lss>2 && !strcmp(ss+2,"LN")))
wcsCelLon_[ii] = jj+1;
else if (!strcmp(ss,"Dec") ||
!strcmp(ss,"b") ||
!strcmp(ss,"Beta") ||
(lss>1 && !strcmp(ss+1,"LAT")) ||
(lss>2 && !strcmp(ss+2,"LT")))
wcsCelLat_[ii] = jj+1;
}
}
// cerr << "wcsCelLon(" << ii << "): " << wcsCelLon_[ii] << endl;
// cerr << "wcsCelLat(" << ii << "): " << wcsCelLat_[ii] << endl;
}
}
}
astEnd;
}
void FitsImage::wcsPhyInit()
{
// now see if we have a 'physical' in WCSP, if so, set LTMV keywords
keyLTMV =0;
char* wcsname = image_->getString("WCSNAMEP");
if (wcsname && *wcsname && !strncmp(wcsname, "PHYSICAL", 8)) {
Vector ll2 = center() - Vector(10,10);
Vector ur2 = center() + Vector(10,10);
double ll[4];
ll[0] =ll2[0];
ll[1] =ll2[1];
ll[2] =1;
ll[3] =1;
double ur[4];
ur[0] =ur2[0];
ur[1] =ur2[1];
ur[2] =1;
ur[3] =1;
setWCSSystem(Coord::WCSP);
int naxes1 = astGetI(astGetFrame(ast_,AST__BASE),"Naxes");
int naxes2 = astGetI(astGetFrame(ast_,AST__CURRENT),"Naxes");
int ss = (naxes1+1)*naxes2;
double* fit = new double[ss];
double tol = 1;
if (astLinearApprox(ast_, ll, ur, tol, fit) != AST__BAD) {
keyLTMV = 1;
imageToPhysical = Matrix(fit[naxes2], fit[naxes2+naxes1],
fit[naxes2+1], fit[naxes2+naxes1+1],
fit[0], fit[1]);
physicalToImage = imageToPhysical.invert();
}
if (fit)
delete [] fit;
}
}
int FitsImage::checkWCS(Vector& vv)
{
// check for reasonable values
return (fabs(vv[0]) < FLT_MAX &&
fabs(vv[1]) < FLT_MAX) ? 1 : 0;
}
int FitsImage::checkWCS(Vector3d& vv)
{
// check for reasonable values
return (fabs(vv[0]) < FLT_MAX &&
fabs(vv[1]) < FLT_MAX &&
fabs(vv[2]) < FLT_MAX ) ? 1 : 0;
}
void FitsImage::setWCSSystem(Coord::CoordSystem sys)
{
if (wcsState_->wcsSystem_ != sys)
if (wcsSystem(ast_,sys))
wcsState_->wcsSystem_ = sys;
}
void FitsImage::setWCSSysSkyFrame(Coord::CoordSystem sys, Coord::SkyFrame sky)
{
if (wcsState_->wcsSystem_ != sys)
if (wcsSystem(ast_,sys))
wcsState_->wcsSystem_ = sys;
if (wcsState_->wcsSkyFrame_ != sky)
if (hasWCSEqu(sys)) {
wcsSkyFrame(ast_,sky);
wcsState_->wcsSkyFrame_ = sky;
}
}
void FitsImage::setWCSFormat(Coord::CoordSystem sys, Coord::SkyFrame sky,
Coord::SkyFormat format)
{
int id = sys-Coord::WCS;
// spacial axes
if (wcsCelLon_[id] && wcsCelLat_[id]) {
switch (format) {
case Coord::DEGREES:
{
ostringstream str;
str << "d." << context_->parent_->precLenDeg_;
wcsFormat(ast_, wcsCelLon_[id], str.str().c_str());
wcsFormat(ast_, wcsCelLat_[id], str.str().c_str());
}
break;
case Coord::SEXAGESIMAL:
{
ostringstream hms;
ostringstream dms;
hms << "hms." << context_->parent_->precHMS_;
dms << "+dms." << context_->parent_->precDMS_;
switch (sky) {
case Coord::FK4:
case Coord::FK5:
case Coord::ICRS:
wcsFormat(ast_, wcsCelLon_[id], hms.str().c_str());
wcsFormat(ast_, wcsCelLat_[id], dms.str().c_str());
break;
case Coord::GALACTIC:
case Coord::ECLIPTIC:
wcsFormat(ast_, wcsCelLon_[id], dms.str().c_str());
wcsFormat(ast_, wcsCelLat_[id], dms.str().c_str());
break;
}
}
break;
}
}
// any other axes
ostringstream str;
str << "%%1." << context_->parent_->precLinear_ << 'G';
for (int ii=0; ii<wcsNaxes_[id]; ii++)
if (!(wcsCelLon_[id] && wcsCelLat_[id]))
wcsFormat(ast_, ii+1, str.str().c_str());
}
Vector FitsImage::vRadToDeg(const Vector& vv, Coord::CoordSystem sys)
{
Vector out = vv;
if (hasWCSCel(sys)) {
int id = sys-Coord::WCS;
for (int ii=0; ii<2; ii++) {
if (wcsCelLon_[id] == ii+1 || wcsCelLat_[id] == ii+1)
out[ii] = radToDeg(out[ii]);
}
}
return out;
}
Vector3d FitsImage::vRadToDeg(const Vector3d& vv, Coord::CoordSystem sys)
{
Vector3d out = vv;
if (hasWCSCel(sys)) {
int id = sys-Coord::WCS;
for (int ii=0; ii<3; ii++) {
if (wcsCelLon_[id] == ii+1 || wcsCelLat_[id] == ii+1)
out[ii] = radToDeg(out[ii]);
}
}
return out;
}
Vector FitsImage::vDegToRad(const Vector& vv, Coord::CoordSystem sys)
{
Vector out = vv;
if (hasWCSCel(sys)) {
int id = sys-Coord::WCS;
for (int ii=0; ii<2; ii++) {
if (wcsCelLon_[id] == ii+1 || wcsCelLat_[id] == ii+1)
out[ii] = degToRad(out[ii]);
}
}
return out;
}
Vector3d FitsImage::vDegToRad(const Vector3d& vv, Coord::CoordSystem sys)
{
Vector3d out = vv;
if (hasWCSCel(sys)) {
int id = sys-Coord::WCS;
for (int ii=0; ii<3; ii++) {
if (wcsCelLon_[id] == ii+1 || wcsCelLat_[id] == ii+1)
out[ii] = degToRad(out[ii]);
}
}
return out;
}
static void fits2TAB(AstFitsChan* chan, const char* extname,
int extver, int extlevel, int* status)
{
FitsFile* ext = ((FitsImage*)astChannelData)->fitsFile();
// just in case
if (!ext) {
*status = 0;
return;
}
// skip the current HDU
ext = new FitsMosaicNextMMapIncr(ext);
while (1) {
// EOF?
if (!ext || !ext->isValid()) {
if (ext)
delete ext;
*status = 0;
return;
}
// is it a bin table?
if (ext->isBinTable()) {
// check its extname
const char* name = ext->extname();
int ver = ext->extver();
int level = ext->extlevel();
if (name && !strcmp(extname,name) && extver==ver && extlevel==level)
break;
}
FitsFile* ptr = ext;
ext = new FitsMosaicNextMMapIncr(ptr);
delete ptr;
}
// ok, found it
astClearStatus; // just to make sure
astBegin; // start memory management
FitsHead* hd = ext->head();
FitsBinTableHDU* hdu = (FitsBinTableHDU*)hd->hdu();
int cols = hdu->cols();
int rows = hdu->rows();
int rowlen = hdu->width();
// create fitstable
AstFitsChan* header = astFitsChan(NULL, NULL, " ");
char* cards = hd->cards();
int ncards = hd->ncard();
for (int ii=0; ii<ncards; ii++) {
char buf[81];
strncpy(buf,cards+(ii*80),80);
buf[80] = '\0';
astPutFits(header, buf, 0);
}
AstFitsTable* table = (AstFitsTable*)astFitsTable(header," ");
for (int ii=0; ii<cols; ii++) {
FitsBinColumn* col = (FitsBinColumn*)hdu->find(ii);
int width = col->width();
int repeat = col->repeat();
char* ptr = (char*)ext->data();
unsigned char* data = new unsigned char[width*rows];
memset(data,0,width*rows);
switch (col->type()) {
case 'I':
for (int ii=0; ii<rows; ii++, ptr+=rowlen)
for (int jj=0; jj<repeat; jj++) {
short vv = col->value(ptr,jj);
memcpy(data+ii*width+jj*2,&vv,2);
}
break;
case 'J':
for (int ii=0; ii<rows; ii++, ptr+=rowlen)
for (int jj=0; jj<repeat; jj++) {
int vv = col->value(ptr,jj);
memcpy(data+ii*width+jj*4,&vv,4);
}
break;
case 'E':
for (int ii=0; ii<rows; ii++, ptr+=rowlen)
for (int jj=0; jj<repeat; jj++) {
float vv = col->value(ptr,jj);
memcpy(data+ii*width+jj*4,&vv,4);
}
break;
case 'D':
for (int ii=0; ii<rows; ii++, ptr+=rowlen)
for (int jj=0; jj<repeat; jj++) {
double vv = col->value(ptr,jj);
memcpy(data+ii*width+jj*8,&vv,8);
}
break;
}
astPutColumnData(table, col->ttype(), 0, width*rows, data);
if (data)
delete [] data;
}
astPutTable(chan, table, extname);
astEnd; // now, clean up memory
if (ext)
delete ext;
*status = 1;
}
AstFrameSet* FitsImage::fits2ast(FitsHead* hd)
{
// we may have an error, just reset
astClearStatus;
astBegin;
// new fitschan
AstFitsChan* chan = astFitsChan(NULL, NULL, " ");
if (!astOK || chan == AST__NULL)
return NULL;
// enable -TAB
astSetI(chan,"TabOK",1);
astPutChannelData(chan, this);
astTableSource(chan, fits2TAB);
// no warning messages
astClear(chan,"Warnings");
// fill up chan
char* cards =NULL;
int ncards =0;
if (hd) {
cards = hd->cards();
ncards = hd->ncard();
}
if (cards == NULL || ncards == 0)
return NULL;
for (int i=0; i<ncards; i++) {
char buf[81];
strncpy(buf,cards+(i*80),80);
buf[80] = '\0';
astPutFits(chan, buf, 0);
// sometimes, we get a bad parse, just ignore
if (!astOK)
astClearStatus;
}
// encoding
// must come after astPutsFits and before astRead
const char* encode = astGetC(chan, "Encoding");
if (encode)
encoding_ = dupstr(encode);
// astSet (chan, "Encoding=FITS-WCS");
// we may have an error, just reset
astClearStatus;
astClear(chan, "Card");
// parse header
AstFrameSet* frameSet = (AstFrameSet*)astRead(chan);
// do we have anything?
if (!astOK || frameSet == AST__NULL ||
strncmp(astGetC(frameSet,"Class"), "FrameSet", 8))
return NULL;
// warn if no inverse
wcsInv_ = astGetI(frameSet, "TranInverse");
if (!wcsInv_)
internalError("Warning: the WCS has no defined inverse. Some functionality may not be available.");
astExport(frameSet);
astEnd;
return frameSet;
}