QuantitiesByAst.hxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// SVN: $Id: QuantitiesByAst.hxx 4813 2008-05-29 22:55:39Z RossCollins $
//-----------------------------------------------------------------------------
#ifndef QUANTITIES_HXX
#define QUANTITIES_HXX

// We will make use of AST for operations involving WCS
extern "C"
{
  #include "ast.h"
}

#include <cmath>
#include <string>

#include "Constants.h"
#include "Converter.h"
#include "FitsFile.h"
#include "Logger.h"
#include "DataMethod.hxx"
#include "DataOps.h"
#include "StringOps.h"
//-----------------------------------------------------------------------------
template<typename DataType>
class QuantitiesByAst : public DataMethod<DataType>
{
public:
  QuantitiesByAst(const FitsFile& aFile)
    : mFile(aFile)
    { DataMethod<DataType>::methodName = "QuantitiesByAst"; }

  virtual ~QuantitiesByAst() { }

  void setSource(IntMap& tbls, IntMap& cols, StringMap units)
    { extnumColno = cols[0]; }

  void setTarget(IntMap& cols, StringMap units)
  {
    raColno     = cols[0];
    decColno    = cols[1];
    posangColno = cols[2];
    extentColno = cols[3];
    xsizeColno  = cols[4];
    ysizeColno  = cols[5];
    minRaColno  = cols[6];
    minDecColno = cols[7];
    maxRaColno  = cols[8];
    maxDecColno = cols[9];

    // AST will calculate angles in radians, so we need the conversion
    // factors for the required units in the schema
    Converter conv;
    raFact     = conv.getFactor("RADIANS", units[0]);
    decFact    = conv.getFactor("RADIANS", units[1]);
    posangFact = conv.getFactor("RADIANS", units[2]);
    extentFact = conv.getFactor("RADIANS", units[3]);
    xsizeFact  = conv.getFactor("RADIANS", units[4]);
    ysizeFact  = conv.getFactor("RADIANS", units[5]);
    minRaFact  = conv.getFactor("RADIANS", units[6]);
    minDecFact = conv.getFactor("RADIANS", units[7]);
    maxRaFact  = conv.getFactor("RADIANS", units[8]);
    maxDecFact = conv.getFactor("RADIANS", units[9]);
  }

  void doit(TableData<DataType>& data, int row1, int row2)
  {
    std::string filename = mFile.getFileName();
    double crval1[4],crval2[4];
    int npixx[4], npixy[4];

    // Use these arrays for AST input
    double xpix[5], ypix[5], xsky[5], ysky[5];

    // Set up AST
    astBegin;

    // Create a new AST fits channel
    AstFitsChan *fitschan = astFitsChan( NULL, NULL, "" );
    AstFrameSet *wcsinfo0 = 0;
    AstFrameSet *wcsinfo1 = 0;
    AstFrameSet *wcsinfo2 = 0;
    AstFrameSet *wcsinfo3 = 0;

    // Get AST data from file
    for (int hduNo = 2; hduNo <= mFile.getNumHdus(); ++hduNo)
    {
      int ext = hduNo - 2;
      // Move to this extension on the FITS file
      mFile.movabsHdu(hduNo);

      // Put all cards in the header of this extension into the AST
      // fits channel
      int ncards = mFile.getNumCards();
      for (int n = 0; n < ncards; ++n) {
          char card[80];
          mFile.readCard(n, card);
          astPutFits(fitschan, card, 0);
      }
      // Make AST wcsinfo object
      astClear( fitschan, "Card" );
      mFile.readKey(FitsFile::cCRVAL1,crval1[ext]);
      mFile.readKey(FitsFile::cCRVAL2,crval2[ext]);
      if (crval1[ext] == 0 && crval2[ext] == 0) {
          log.addError("Could not read WCS info for " + filename);
      } else {
        if (ext == 0) {
          wcsinfo0 = (AstFrameSet*)astRead(fitschan);
        } else if (ext == 1) {
          wcsinfo1 = (AstFrameSet*)astRead(fitschan);
        } else if (ext == 2) {
          wcsinfo2 = (AstFrameSet*)astRead(fitschan);
        } else if (ext == 3) {
          wcsinfo3 = (AstFrameSet*)astRead(fitschan);
        }
        // Centre of image in pixel coordinates

        if (mFile.hasKey("ZIMAGE"))
        {
          mFile.readKey(FitsFile::ZNAXIS1, npixx[ext]);
          mFile.readKey(FitsFile::ZNAXIS2, npixy[ext]);
        }
        else
        {
          mFile.readKey(FitsFile::NAXIS1, npixx[ext]);
          mFile.readKey(FitsFile::NAXIS2, npixy[ext]);
        }
      }
    }

    for (int rowNo = row1; rowNo < row2; ++rowNo)
    {
      // Get the FITS extension number for the data in this row. Note
      // that the extension number is assumed to be in the fitsio
      // convention
      Numeric value;
      data.value(extnumColno, rowNo, value);
      int extNum = (int)value - 2;

      if (crval1[extNum] != 0 || crval2[extNum] != 0)
      {
        // Not sure how pixel coordinates are defined. What is (x,y) of the
        // centre of the first pixel?
        xpix[0] = (double)npixx[extNum] / 2;
        ypix[0] = (double)npixy[extNum] / 2;

        // Pixel coordinates of the four corners
        xpix[1] = 0;
        ypix[1] = 0;
        xpix[2] = npixx[extNum];
        ypix[2] = 0;
        xpix[3] = 0;
        ypix[3] = npixy[extNum];
        xpix[4] = npixx[extNum];
        ypix[4] = npixy[extNum];

        // Get the sky coordinates for the array of pixel coordinates. Note
        // these will be in radians. Angles can also be negative
        if (extNum == 0) {
          astTran2(wcsinfo0, 5, xpix, ypix, 1, xsky, ysky);
        } else if (extNum == 1) {
          astTran2(wcsinfo1, 5, xpix, ypix, 1, xsky, ysky);
        } else if (extNum == 2) {
          astTran2(wcsinfo2, 5, xpix, ypix, 1, xsky, ysky);
        } else if (extNum == 3) {
          astTran2(wcsinfo3, 5, xpix, ypix, 1, xsky, ysky);
        }

        // Make sure all RAs are in the range 0 to 2*pi
        for (int n = 0; n < 5; ++n)
        {
          if (xsky[n] < 0)
          {
            xsky[n] += 2*M_PI;
          }
        }
        // Get the central ra and dec in units as required in the schema
        Numeric raCen  = raFact * xsky[0];
        Numeric decCen = decFact * ysky[0];

        // Calculate max and min of RA and DEC from xsky, ysky
        // first find max and min values
        double xskyMin = 9999.999;
        double xskyMax = -9999.999;
        double yskyMin = 9999.999;
        double yskyMax = -9999.999;
        for (int n = 1; n < 5; ++n)
        {
          if (xsky[n] < xskyMin) xskyMin = xsky[n];
          if (xsky[n] > xskyMax) xskyMax = xsky[n];
          if (ysky[n] < yskyMin) yskyMin = ysky[n];
          if (ysky[n] > yskyMax) yskyMax = ysky[n];
        }
        // then project RA in 0 to 3*pi for fields crossing RA=0
        double h = 0.0;
        if (xskyMax - xskyMin > M_PI)
        {
          h = xskyMax;
          xskyMax = xskyMin + 2*M_PI;
          xskyMin = h;
        }
        Numeric raMin  = minRaFact * xskyMin;
        Numeric decMin = minDecFact * yskyMin;
        Numeric raMax  = maxRaFact * xskyMax;
        Numeric decMax = maxDecFact * yskyMax;

        // Find angular separation between the centre and each of the
        // four corners and find the maximum of these
        double maxsep = 0;
        for (int n = 0; n < 5; ++n)
        {
          double ang = DataOps::calcAngsep(xsky[0], ysky[0],
                                            xsky[n], ysky[n]);
          if (ang > maxsep) maxsep = ang;
        }
        // Take image extent as maximum separation and put in required units
        Numeric imageExtent = maxsep * extentFact;

        // Find the posAngle
        // Define N-S line by incrementing the ysky coordinate (usually dec)
        // by a small amount
        xsky[1] = xsky[0];
        ysky[1] = ysky[0] + 0.01;
        // Find equiavelent pixel corrinates for the N-S point
        if (extNum == 0) {
          astTran2(wcsinfo0, 5, xsky, ysky, 0, xpix, ypix);
        } else if (extNum == 1) {
          astTran2(wcsinfo1, 5, xsky, ysky, 0, xpix, ypix);
        } else if (extNum == 2) {
          astTran2(wcsinfo2, 5, xsky, ysky, 0, xpix, ypix);
        } else if (extNum == 3) {
          astTran2(wcsinfo3, 5, xsky, ysky, 0, xpix, ypix);
        }
        // Determine angle between N-S line and positive x-axis with origin
        // at the centre of the image
        Numeric posAng = posangFact * atan2(ypix[1] - ypix[0],
                                            xpix[1] - xpix[0]);

        // X and Y pixel size is calculated by incrementing the
        // respective pixel sizes at the image centre
        xpix[1] = xpix[0] + 1;
        ypix[1] = ypix[0];
        xpix[2] = xpix[0];
        ypix[2] = ypix[0] + 1;
        if (extNum == 0) {
          astTran2(wcsinfo0, 5, xpix, ypix, 1, xsky, ysky);
        } else if (extNum == 1) {
          astTran2(wcsinfo1, 5, xpix, ypix, 1, xsky, ysky);
        } else if (extNum == 2) {
          astTran2(wcsinfo2, 5, xpix, ypix, 1, xsky, ysky);
        } else if (extNum == 3) {
          astTran2(wcsinfo3, 5, xpix, ypix, 1, xsky, ysky);
        }
        Numeric xPixSize = xsizeFact * DataOps::calcAngsep(xsky[0], ysky[0],
                                                           xsky[1], ysky[1]);
        Numeric yPixSize = ysizeFact * DataOps::calcAngsep(xsky[0], ysky[0],
                                                           xsky[2], ysky[2]);

        data.assign(raColno, rowNo, raCen);
        data.assign(decColno, rowNo, decCen);
        data.assign(minRaColno, rowNo, raMin);
        data.assign(minDecColno, rowNo, decMin);
        data.assign(maxRaColno, rowNo, raMax);
        data.assign(maxDecColno, rowNo, decMax);
        data.assign(extentColno, rowNo, imageExtent);
        data.assign(posangColno, rowNo, posAng);
        data.assign(xsizeColno, rowNo, xPixSize);
        data.assign(ysizeColno, rowNo, yPixSize);
      }
    }
    // Finshed with AST
    astEnd;
  }
  // No data from raw data included
  void doit2(TableData<DataType>& rawData, TableData<DataType>& data, int row1, int row2) { }

private:
  FitsFile mFile;
  Logger log;

  int extnumColno, raColno, decColno, posangColno;
  int extentColno, xsizeColno, ysizeColno;
  int minRaColno,maxRaColno,minDecColno,maxDecColno;
  double raFact, decFact, posangFact, extentFact, xsizeFact, ysizeFact;
  double minRaFact, maxRaFact, minDecFact, maxDecFact;
};
//-----------------------------------------------------------------------------
#endif
//-----------------------------------------------------------------------------
// Change log:
//
//  7-Jun-2005, ETWS: Original version.
// 22-Mar-2007, RSC:  Robustified against dodgy WCS values.
//  2-May-2007,  RSC: Reverted to original design of allowing dodgy WCS values.
//  7-Apr-2008, ETWS: Upgraded to use new detection table layout.

---