PeriodicBinaryTransform.cxx

Go to the documentation of this file.

// for wcslib
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include "transforms/PeriodicBinaryTransform.h"

#include "axes/AxisModelBase.h"
#include "graphics/Rectangle.h"
#include "UnaryTransform.h"

#ifdef HAVE_WCSLIB
#include <string.h>
#ifdef _WIN32
#include "C/wcs.h"
#else
#include "wcslib/wcs.h"
#endif
#endif

#include <cmath>
#include <cstdio>
#include <stdexcept>
#include <sstream>
#include <cassert>

using namespace hippodraw;

using std::max;
using std::abs;
using std::vector;

PeriodicBinaryTransform::PeriodicBinaryTransform() :
  BinaryTransform(),
  m_x_limits ( -180.0, +180.0 ),
  m_y_limits ( - 90.0, + 90.0 ),
  m_x_offset ( 0.0 ),
  m_y_offset ( 0.0 )
{
}

PeriodicBinaryTransform::PeriodicBinaryTransform ( UnaryTransform * z,
                                                   bool is_periodic,
                                                   bool needs_grid,
                                                   bool needs_x_ticks,
                                                   bool needs_y_ticks,
                                                   double xlo, double xhi,
                                                   double ylo, double yhi) :
  BinaryTransform ( z, is_periodic, needs_grid, needs_x_ticks, needs_y_ticks ),
  m_x_limits ( xlo, xhi ),
  m_y_limits ( ylo, yhi ),
  m_x_offset ( 0.0 ),
  m_y_offset ( 0.0 )
{
}

PeriodicBinaryTransform::
PeriodicBinaryTransform ( const PeriodicBinaryTransform & t ) :
  BinaryTransform ( t ),
  m_x_limits( t.limitX() ),
  m_y_limits( t.limitY() ),
  m_x_offset( t.xOffset() ),
  m_y_offset( t.yOffset() )
{
}

PeriodicBinaryTransform::~PeriodicBinaryTransform ()
{
  // Does nothing 
}

const Range & PeriodicBinaryTransform::limitX () const
{
  return m_x_limits;
}


const Range & PeriodicBinaryTransform::limitY () const
{
  return m_y_limits;
}


double PeriodicBinaryTransform::xOffset() const
{
  return m_x_offset;
}

void PeriodicBinaryTransform::setXOffset( double x_offset ) 
{
  m_x_offset = x_offset;
}

double PeriodicBinaryTransform::yOffset() const
{
  return m_y_offset;
}

void PeriodicBinaryTransform::setYOffset( double y_offset )
{
  
  m_y_offset = y_offset;
}


double
PeriodicBinaryTransform::
moduloAdd( double a1, double a2, hippodraw::Axes::Type axis ) const
{
  if ( axis == Axes::X ) {
    return moduloAddX ( a1,  a2);
  }
  else if ( axis == Axes::Y ) {
    return moduloAddY ( a1,  a2);
  }
  else  return a1 + a2;
}

double
PeriodicBinaryTransform::
moduloSub( double s1, double s2, hippodraw::Axes::Type axis ) const
{
  if ( axis == Axes::X ) {
    return moduloSubX ( s1,  s2);
  }
  else if ( axis == Axes::Y ) {
    return moduloSubY ( s1,  s2);
  }
  else return s1 - s2;
}

double PeriodicBinaryTransform::moduloAddX( double x1, double x2 ) const
{
  if ( x2 < -DBL_EPSILON )
    return moduloSubX ( x1,  -x2 );
  
  double overshoot = x1 + x2 - m_x_limits.high();
  
  if ( overshoot > DBL_EPSILON ) {
    return  m_x_limits.low() + overshoot;
  }
  else {
    return x1 + x2;  
  }
}

double PeriodicBinaryTransform::moduloSubX( double x1, double x2 ) const
{
  if ( x2 < -DBL_EPSILON )
    return moduloAddX ( x1, -x2 );
  
  double undershoot = m_x_limits.low() - ( x1 - x2 );

  if ( undershoot > DBL_EPSILON) 
    return  m_x_limits.high() - undershoot;
  else
    return x1 - x2;      
}


double PeriodicBinaryTransform::moduloAddY( double y1, double y2 ) const
{
  if ( y2 < -DBL_EPSILON ) return moduloSubY ( y1,  -y2 );
  
  double overshoot = y1 + y2 - m_y_limits.high();
    
  if ( overshoot > DBL_EPSILON ) {
    return  m_y_limits.low() + overshoot;
  }
  else {
    return y1 + y2;
  }
  
}

double PeriodicBinaryTransform::moduloSubY( double y1, double y2 ) const
{
  if ( y2 < -DBL_EPSILON ) return moduloAddY ( y1, -y2 );
  
  double undershoot = m_y_limits.low() - ( y1 - y2 );
    
  if ( undershoot > DBL_EPSILON )  return  m_y_limits.high() - undershoot;
  else return y1 - y2;
}


Rect PeriodicBinaryTransform::calcRectangle ( const Range & lat,
                                           const Range & lon )
{
  double x_lo = lat.low ();
  double x_hi = lat.high ();
  
  double y_lo = lon.low ();
  double y_hi = lon.high ();

  double x, y;
  
  double x_min =  1000;
  double x_max = -1000;
  double y_min =  1000;
  double y_max = -1000;

  int n = 50;
  double dx = ( x_hi - x_lo ) / n;
  double dy = ( y_hi - y_lo ) / n;

  
  // Finding out xmin, xmax, ymin, ymax along line y =  y_lo
  for ( int i = 0; i <= n; i++)
    {
      x = x_lo + i * dx; 
      y = y_lo;
      transform ( x, y );
      x_min = ( x_min < x ) ? x_min : x;
      x_max = ( x_max > x ) ? x_max : x;
      y_min = ( y_min < y ) ? y_min : y;
      y_max = ( y_max > y ) ? y_max : y;
    }

  // Finding out xmin, xmax, ymin, ymax along line y =  y_hi
  for ( int i = 0; i <= n; i++)
    {
      x = x_lo + i * dx; 
      y = y_hi;
      transform ( x, y );
      x_min = ( x_min < x ) ? x_min : x;
      x_max = ( x_max > x ) ? x_max : x;
      y_min = ( y_min < y ) ? y_min : y;
      y_max = ( y_max > y ) ? y_max : y;
    }

  // Finding out xmin, xmax, ymin, ymax along line x =  x_lo
  for ( int i = 0; i <= n; i++)
    {
      x = x_lo; 
      y = y_lo + i * dy;
      transform ( x, y );
      x_min = ( x_min < x ) ? x_min : x;
      x_max = ( x_max > x ) ? x_max : x;
      y_min = ( y_min < y ) ? y_min : y;
      y_max = ( y_max > y ) ? y_max : y;
    }

  // Finding out xmin, xmax, ymin, ymax along line x =  x_hi
  for ( int i = 0; i <= n; i++)
    {
      x = x_hi; 
      y = y_lo + i * dy;
      transform ( x, y );
      x_min = ( x_min < x ) ? x_min : x;
      x_max = ( x_max > x ) ? x_max : x;
      y_min = ( y_min < y ) ? y_min : y;
      y_max = ( y_max > y ) ? y_max : y;
    }
    
  return Rect ( x_min, y_min, x_max - x_min, y_max - y_min );
}


void PeriodicBinaryTransform::validate ( Range & lat, Range & lon ) const
{
  if ( lat.low () < m_x_limits.low() ) lat.setLow ( m_x_limits.low() );
  if ( lat.high () > m_x_limits.high() ) lat.setHigh ( m_x_limits.high() );

  if ( lon.low () < m_y_limits.low() ) lon.setLow ( m_y_limits.low() );
  if ( lon.high () > m_y_limits.high() ) lon.setHigh ( m_y_limits.high() );
}

const vector < AxisTick > &
PeriodicBinaryTransform::
setTicks ( AxisModelBase & model, hippodraw::Axes::Type axis )
{
  if ( axis == Axes::Z ) {
    return m_z -> setTicks ( model );
  }
  //else
  setTickStep ( model );
  setFirstTick ( model );
  return genTicks ( model, axis );
}

void
PeriodicBinaryTransform::
adjustValues ( AxisModelBase &,
               hippodraw::Axes::Type,
               const Range & )
{
  // Does not do anything as of now 
  return;
}

inline double FLT_EQUAL( double x, double y )
{
  return ( (double)abs( x - y ) <= 2.0 * ( y * FLT_EPSILON + FLT_MIN ) );
}


void PeriodicBinaryTransform::setTickStep( AxisModelBase & axis )
{
  const Range & range = axis.getRange(false);
  double rangeLength = range.length();

  double scale_factor = axis.getScaleFactor();
  rangeLength *= scale_factor;
  
  // The following algorithm determines the magnitude of the range...
  double rmag = floor( log10( rangeLength ) );
  axis.setRMag( rmag );
  
  double scalelow = range.low() * scale_factor;
  double scalehigh = range.high() * scale_factor;

  // We will also need the largest magnitude for labels.
  double pmag = max( floor( log10( abs ( scalehigh ) ) ), 
                     floor( log10( abs ( scalelow ) ) ) );
  axis.setPMag( pmag );
  
  axis.setTickStep( rangeLength / 4.0 );
}

void PeriodicBinaryTransform::setFirstTick( AxisModelBase & axis )
{
  const Range & range = axis.getRange(false);
  
  axis.setFirstTick ( range.low() );
}


const vector < AxisTick > &
PeriodicBinaryTransform::
genTicks( AxisModelBase & axis, hippodraw::Axes::Type axistype )
{ 
  double y = 0.0, ylabel;
  
  int num_ticks = 0;
  m_ticks.clear();
  double pmag = axis.getPMag();
  double rmag = axis.getRMag();
  double first_tick = axis.getFirstTick();
  double tick_step  = axis.getTickStep();
  double scale_factor = axis.getScaleFactor();
  
  // pmag will get set to 0 if it is less than or equal to 3.  This
  // is used later to determine scientific notation.  However, m_rmag
  // is still needed as the original magnitude for calculations such
  // as decimal place notation, and rounding to nice numbers.
  
  bool use_pmag = abs ( pmag ) > 3.0;

  axis.setUsePMag ( use_pmag );
  
  char pstr[10];
  char labl[10];

  int decimals = 0;

  // The following if-else block decides the pstr string, which holds
  // the number of decimal places in the label.

  //   if( fabs( m_pmag ) > 3.0 ) {
  if ( use_pmag ) {  
    // If we are greater than mag 3, we are showing scientific
    // notation.  How many decimals we show is determined by the
    // difference between the range magnitude and the power magnitude.
  
    decimals = static_cast<int>( pmag - rmag );
    // minumum 1 decimal in scientific notation
    
    if( !decimals ) decimals++;
  
  } else {
    
    if( rmag > 0.0 ){
    
      // If we are less than mag 3 and positive, then no decimal
      // accuracy is needed.
      
      decimals = 0;
      
    } else {
    
      // If we are less than mag 3 and negative, then we are suddenly
      // looking at decimal numbers not in scientific notation.
      // Therefore we hold as many decimal places as the magnitude.
      
      decimals = static_cast<int>( abs( rmag ) );
      
    }
  
  }
  // @todo decimals should never be negative here, but it does end up
  //    negative in some cases. See the "dirty fix" in Range.cxx, that
  //    dirty-fixed this problem too. But a better fix is needed. 
  if (decimals < 0) 
    decimals = 0;
   
  sprintf( pstr, "%%1.%df", decimals );

  y = first_tick;
  const Range & range = axis.getRange(false);
  double range_high = range.high();
  range_high *= scale_factor;
  
  while( y <= range_high || FLT_EQUAL( range_high, y ) )
    {
      double value = 0.0;
      
      if ( axistype == Axes::X ) {
        value = moduloAddX( y, xOffset() );
      }
      else if ( axistype == Axes::Y ) {
        value = moduloAddY( y, yOffset() );
      }
      
      // Now we either keep the original magnitude
      // or reduce it in order to express it in scientific notation.
      
      if ( use_pmag ) ylabel = value / pow( 10.0, pmag );
      else ylabel = value;
      
      value /= scale_factor;
      sprintf( labl, pstr, ylabel );
      m_ticks.push_back( AxisTick ( y, labl ) );
      
      num_ticks++;
      if ( tick_step == 0.0 ) break;
      y += tick_step;
      
    }
  
  return m_ticks;
}


bool
PeriodicBinaryTransform::
isLinearInXY () const
{
  return false;
}





/* Use WCSLIB to do transform. */

void 
PeriodicBinaryTransform::
initwcs(const std::string &transformName, double* crpix, 
        double* crval, double* cdelt, double crota2, bool galactic)
{
#ifdef HAVE_WCSLIB
    // Assert the sizeof wcsprm struct is smaller than the size allocated. 
    assert(sizeof(wcsprm)<=2000);
    m_wcs = reinterpret_cast<wcsprm*>(m_wcs_struct);
    m_wcs->flag = -1;

    // Call wcsini() in WCSLIB.
    int naxis = 2;
    wcsini(1, naxis, m_wcs);

    // Transfer parameters from c++ class to c stuct.
    std::string 
        lon_type = (galactic? "GLON-" : "RA---") + transformName,
        lat_type =  (galactic? "GLAT-" : "DEC--") + transformName;
    strcpy(m_wcs->ctype[0], lon_type.c_str() );
    strcpy(m_wcs->ctype[1], lat_type.c_str() );

    for( int i=0; i<naxis; ++i){
        m_wcs->crval[i] = crval[i]; // reference value
        m_wcs->crpix[i] = crpix[i]; // pixel coordinate
        m_wcs->cdelt[i] = cdelt[i]; // scale factor
    }
#else
    throwWCSMissing ();
#endif

}

void
PeriodicBinaryTransform::
throwWCSMissing () const
{
  const std::string what
    ( "HippoDraw has not been compiled with WCSLIB support" );
  throw std::runtime_error ( what );
}



/* virtual */
void
PeriodicBinaryTransform::
transform ( double & lon, double & lat ) const
{

#ifdef HAVE_WCSLIB
    int ncoords = 1;
    int nelem = 2;
    double  imgcrd[2], pixcrd[2];
    double phi[1], theta[1];
    int stat[1];

    double worldcrd[] ={lon,lat};

    //    int returncode = 
    wcss2p(m_wcs, ncoords, nelem, worldcrd, phi, theta, imgcrd, pixcrd, stat);
    
    lon = pixcrd[0];
    lat = pixcrd[1]; 
#else
    throwWCSMissing ();
#endif


}


bool
PeriodicBinaryTransform::
inverseTransform ( double & lon, double & lat ) const
{
  
#ifdef HAVE_WCSLIB
    int ncoords = 1;
    int nelem = 2;
    double worldcrd[2], imgcrd[2];
    double phi[1], theta[1];
    int stat[1];

    double pixcrd[] = {lon,lat};

    int returncode = 
    wcsp2s(m_wcs, ncoords, nelem, pixcrd, imgcrd, phi, theta, worldcrd, stat);
    

    // If there is invalid coordinate, return false.
    if ( returncode == 8 ) return false;

    // Make sure the inverse transform result in (-180,180) or (0,360)
    //if (worldcrd[0]<0) worldcrd[0]=worldcrd[0]+360;
    lon = worldcrd [0];
    lat = worldcrd [1];    
        

    return true;
#else
    throwWCSMissing ();
    return false;
#endif
}

void
PeriodicBinaryTransform::
transform ( std::vector< double > & lon,
            std::vector< double > & lat ) const
{     

#ifdef HAVE_WCSLIB
  assert ( lat.size() == lon.size() );
  int size = lat.size();

  // Form vector to array.
  vector < double >  worldcrd ( 2*size );
  for ( int i = 0; i < size; i++ ) {
    worldcrd[2*i]= lon[i];
    worldcrd[2*i+1]= lat[i];
  }
  
  int ncoords = size;
  int nelem = 2;

  vector < double > pixcrd( 2*size );
  vector < double > imgcrd( 2*size );
  vector < double > phi( size );
  vector < double > theta( size );
  vector < int > stat ( size );

  //int returncode =
  wcss2p ( m_wcs, ncoords, nelem, &worldcrd[0], &phi[0], &theta[0],
           &imgcrd[0], &pixcrd[0], &stat[0] );
  
  // From array to vector;  
  for ( int i = 0; i < size; i++ ) {
    lon[i]=pixcrd[2*i];
    lat[i]=pixcrd[2*i+1];
  }
   
#else
    throwWCSMissing ();
#endif
}