PeriodicBinaryTransform.cxx

Go to the documentation of this file.
00001 
00011 // for wcslib
00012 #ifdef HAVE_CONFIG_H
00013 #include "config.h"
00014 #endif
00015 
00016 #include "transforms/PeriodicBinaryTransform.h"
00017 
00018 #include "axes/AxisModelBase.h"
00019 #include "graphics/Rectangle.h"
00020 #include "UnaryTransform.h"
00021 
00022 #ifdef HAVE_WCSLIB
00023 #ifdef _WIN32
00024 #include "wcslib/C/wcs.h"
00025 #else
00026 #include "wcslib/wcs.h"
00027 #endif
00028 #endif
00029 
00030 #include <cmath>
00031 #include <cstdio>
00032 #include <stdexcept>
00033 #include <sstream>
00034 #include <cassert>
00035 
00036 using namespace hippodraw;
00037 
00038 using std::max;
00039 using std::abs;
00040 using std::vector;
00041 
00047 PeriodicBinaryTransform::PeriodicBinaryTransform() :
00048   BinaryTransform(),
00049   m_x_limits ( -180.0, +180.0 ),
00050   m_y_limits ( - 90.0, + 90.0 ),
00051   m_x_offset ( 0.0 ),
00052   m_y_offset ( 0.0 )
00053 {
00054 }
00055 
00056 PeriodicBinaryTransform::PeriodicBinaryTransform ( UnaryTransform * z,
00057                                                    bool is_periodic,
00058                                                    bool needs_grid,
00059                                                    bool needs_x_ticks,
00060                                                    bool needs_y_ticks,
00061                                                    double xlo, double xhi,
00062                                                    double ylo, double yhi) :
00063   BinaryTransform ( z, is_periodic, needs_grid, needs_x_ticks, needs_y_ticks ),
00064   m_x_limits ( xlo, xhi ),
00065   m_y_limits ( ylo, yhi ),
00066   m_x_offset ( 0.0 ),
00067   m_y_offset ( 0.0 )
00068 {
00069 }
00070 
00071 PeriodicBinaryTransform::
00072 PeriodicBinaryTransform ( const PeriodicBinaryTransform & t ) :
00073   BinaryTransform ( t ),
00074   m_x_limits( t.limitX() ),
00075   m_y_limits( t.limitY() ),
00076   m_x_offset( t.xOffset() ),
00077   m_y_offset( t.yOffset() )
00078 {
00079 }
00080 
00081 PeriodicBinaryTransform::~PeriodicBinaryTransform ()
00082 {
00083   // Does nothing 
00084 }
00085 
00086 const Range & PeriodicBinaryTransform::limitX () const
00087 {
00088   return m_x_limits;
00089 }
00090 
00091 
00092 const Range & PeriodicBinaryTransform::limitY () const
00093 {
00094   return m_y_limits;
00095 }
00096 
00097 
00098 double PeriodicBinaryTransform::xOffset() const
00099 {
00100   return m_x_offset;
00101 }
00102 
00103 void PeriodicBinaryTransform::setXOffset( double x_offset ) 
00104 {
00105   m_x_offset = x_offset;
00106 }
00107 
00108 double PeriodicBinaryTransform::yOffset() const
00109 {
00110   return m_y_offset;
00111 }
00112 
00113 void PeriodicBinaryTransform::setYOffset( double y_offset )
00114 {
00115   
00116   m_y_offset = y_offset;
00117 }
00118 
00119 
00120 double
00121 PeriodicBinaryTransform::
00122 moduloAdd( double a1, double a2, hippodraw::Axes::Type axis ) const
00123 {
00124   if ( axis == Axes::X ) {
00125     return moduloAddX ( a1,  a2);
00126   }
00127   else if ( axis == Axes::Y ) {
00128     return moduloAddY ( a1,  a2);
00129   }
00130   else  return a1 + a2;
00131 }
00132 
00133 double
00134 PeriodicBinaryTransform::
00135 moduloSub( double s1, double s2, hippodraw::Axes::Type axis ) const
00136 {
00137   if ( axis == Axes::X ) {
00138     return moduloSubX ( s1,  s2);
00139   }
00140   else if ( axis == Axes::Y ) {
00141     return moduloSubY ( s1,  s2);
00142   }
00143   else return s1 - s2;
00144 }
00145 
00146 double PeriodicBinaryTransform::moduloAddX( double x1, double x2 ) const
00147 {
00148   if ( x2 < -DBL_EPSILON )
00149     return moduloSubX ( x1,  -x2 );
00150   
00151   double overshoot = x1 + x2 - m_x_limits.high();
00152   
00153   if ( overshoot > DBL_EPSILON ) {
00154     return  m_x_limits.low() + overshoot;
00155   }
00156   else {
00157     return x1 + x2;  
00158   }
00159 }
00160 
00161 double PeriodicBinaryTransform::moduloSubX( double x1, double x2 ) const
00162 {
00163   if ( x2 < -DBL_EPSILON )
00164     return moduloAddX ( x1, -x2 );
00165   
00166   double undershoot = m_x_limits.low() - ( x1 - x2 );
00167 
00168   if ( undershoot > DBL_EPSILON) 
00169     return  m_x_limits.high() - undershoot;
00170   else
00171     return x1 - x2;      
00172 }
00173 
00174 
00175 double PeriodicBinaryTransform::moduloAddY( double y1, double y2 ) const
00176 {
00177   if ( y2 < -DBL_EPSILON ) return moduloSubY ( y1,  -y2 );
00178   
00179   double overshoot = y1 + y2 - m_y_limits.high();
00180     
00181   if ( overshoot > DBL_EPSILON ) {
00182     return  m_y_limits.low() + overshoot;
00183   }
00184   else {
00185     return y1 + y2;
00186   }
00187   
00188 }
00189 
00190 double PeriodicBinaryTransform::moduloSubY( double y1, double y2 ) const
00191 {
00192   if ( y2 < -DBL_EPSILON ) return moduloAddY ( y1, -y2 );
00193   
00194   double undershoot = m_y_limits.low() - ( y1 - y2 );
00195     
00196   if ( undershoot > DBL_EPSILON )  return  m_y_limits.high() - undershoot;
00197   else return y1 - y2;
00198 }
00199 
00200 
00201 Rect PeriodicBinaryTransform::calcRectangle ( const Range & lat,
00202                                            const Range & lon )
00203 {
00204   double x_lo = lat.low ();
00205   double x_hi = lat.high ();
00206   
00207   double y_lo = lon.low ();
00208   double y_hi = lon.high ();
00209 
00210   double x, y;
00211   
00212   double x_min =  1000;
00213   double x_max = -1000;
00214   double y_min =  1000;
00215   double y_max = -1000;
00216 
00217   int n = 50;
00218   double dx = ( x_hi - x_lo ) / n;
00219   double dy = ( y_hi - y_lo ) / n;
00220 
00221   
00222   // Finding out xmin, xmax, ymin, ymax along line y =  y_lo
00223   for ( int i = 0; i <= n; i++)
00224     {
00225       x = x_lo + i * dx; 
00226       y = y_lo;
00227       transform ( x, y );
00228       x_min = ( x_min < x ) ? x_min : x;
00229       x_max = ( x_max > x ) ? x_max : x;
00230       y_min = ( y_min < y ) ? y_min : y;
00231       y_max = ( y_max > y ) ? y_max : y;
00232     }
00233 
00234   // Finding out xmin, xmax, ymin, ymax along line y =  y_hi
00235   for ( int i = 0; i <= n; i++)
00236     {
00237       x = x_lo + i * dx; 
00238       y = y_hi;
00239       transform ( x, y );
00240       x_min = ( x_min < x ) ? x_min : x;
00241       x_max = ( x_max > x ) ? x_max : x;
00242       y_min = ( y_min < y ) ? y_min : y;
00243       y_max = ( y_max > y ) ? y_max : y;
00244     }
00245 
00246   // Finding out xmin, xmax, ymin, ymax along line x =  x_lo
00247   for ( int i = 0; i <= n; i++)
00248     {
00249       x = x_lo; 
00250       y = y_lo + i * dy;
00251       transform ( x, y );
00252       x_min = ( x_min < x ) ? x_min : x;
00253       x_max = ( x_max > x ) ? x_max : x;
00254       y_min = ( y_min < y ) ? y_min : y;
00255       y_max = ( y_max > y ) ? y_max : y;
00256     }
00257 
00258   // Finding out xmin, xmax, ymin, ymax along line x =  x_hi
00259   for ( int i = 0; i <= n; i++)
00260     {
00261       x = x_hi; 
00262       y = y_lo + i * dy;
00263       transform ( x, y );
00264       x_min = ( x_min < x ) ? x_min : x;
00265       x_max = ( x_max > x ) ? x_max : x;
00266       y_min = ( y_min < y ) ? y_min : y;
00267       y_max = ( y_max > y ) ? y_max : y;
00268     }
00269     
00270   return Rect ( x_min, y_min, x_max - x_min, y_max - y_min );
00271 }
00272 
00273 
00274 void PeriodicBinaryTransform::validate ( Range & lat, Range & lon ) const
00275 {
00276   if ( lat.low () < m_x_limits.low() ) lat.setLow ( m_x_limits.low() );
00277   if ( lat.high () > m_x_limits.high() ) lat.setHigh ( m_x_limits.high() );
00278 
00279   if ( lon.low () < m_y_limits.low() ) lon.setLow ( m_y_limits.low() );
00280   if ( lon.high () > m_y_limits.high() ) lon.setHigh ( m_y_limits.high() );
00281 }
00282 
00286 const vector < AxisTick > &
00287 PeriodicBinaryTransform::
00288 setTicks ( AxisModelBase & model, hippodraw::Axes::Type axis )
00289 {
00290   if ( axis == Axes::Z ) {
00291     return m_z -> setTicks ( model );
00292   }
00293   //else
00294   setTickStep ( model );
00295   setFirstTick ( model );
00296   return genTicks ( model, axis );
00297 }
00298 
00300 void
00301 PeriodicBinaryTransform::
00302 adjustValues ( AxisModelBase &,
00303                hippodraw::Axes::Type,
00304                const Range & )
00305 {
00306   // Does not do anything as of now 
00307   return;
00308 }
00309 
00310 inline double FLT_EQUAL( double x, double y )
00311 {
00312   return ( (double)abs( x - y ) <= 2.0 * ( y * FLT_EPSILON + FLT_MIN ) );
00313 }
00314 
00315 
00316 void PeriodicBinaryTransform::setTickStep( AxisModelBase & axis )
00317 {
00318   const Range & range = axis.getRange(false);
00319   double rangeLength = range.length();
00320 
00321   double scale_factor = axis.getScaleFactor();
00322   rangeLength *= scale_factor;
00323   
00324   // The following algorithm determines the magnitude of the range...
00325   double rmag = floor( log10( rangeLength ) );
00326   axis.setRMag( rmag );
00327   
00328   double scalelow = range.low() * scale_factor;
00329   double scalehigh = range.high() * scale_factor;
00330 
00331   // We will also need the largest magnitude for labels.
00332   double pmag = max( floor( log10( abs ( scalehigh ) ) ), 
00333                      floor( log10( abs ( scalelow ) ) ) );
00334   axis.setPMag( pmag );
00335   
00336   axis.setTickStep( rangeLength / 4.0 );
00337 }
00338 
00339 void PeriodicBinaryTransform::setFirstTick( AxisModelBase & axis )
00340 {
00341   const Range & range = axis.getRange(false);
00342   
00343   axis.setFirstTick ( range.low() );
00344 }
00345 
00346 
00351 const vector < AxisTick > &
00352 PeriodicBinaryTransform::
00353 genTicks( AxisModelBase & axis, hippodraw::Axes::Type axistype )
00354 { 
00355   double y = 0.0, ylabel;
00356   
00357   int num_ticks = 0;
00358   m_ticks.clear();
00359   double pmag = axis.getPMag();
00360   double rmag = axis.getRMag();
00361   double first_tick = axis.getFirstTick();
00362   double tick_step  = axis.getTickStep();
00363   double scale_factor = axis.getScaleFactor();
00364   
00365   // pmag will get set to 0 if it is less than or equal to 3.  This
00366   // is used later to determine scientific notation.  However, m_rmag
00367   // is still needed as the original magnitude for calculations such
00368   // as decimal place notation, and rounding to nice numbers.
00369   
00370   bool use_pmag = abs ( pmag ) > 3.0;
00371 
00372   axis.setUsePMag ( use_pmag );
00373   
00374   char pstr[10];
00375   char labl[10];
00376 
00377   int decimals = 0;
00378 
00379   // The following if-else block decides the pstr string, which holds
00380   // the number of decimal places in the label.
00381 
00382   //   if( fabs( m_pmag ) > 3.0 ) {
00383   if ( use_pmag ) {  
00384     // If we are greater than mag 3, we are showing scientific
00385     // notation.  How many decimals we show is determined by the
00386     // difference between the range magnitude and the power magnitude.
00387   
00388     decimals = static_cast<int>( pmag - rmag );
00389     // minumum 1 decimal in scientific notation
00390     
00391     if( !decimals ) decimals++;
00392   
00393   } else {
00394     
00395     if( rmag > 0.0 ){
00396     
00397       // If we are less than mag 3 and positive, then no decimal
00398       // accuracy is needed.
00399       
00400       decimals = 0;
00401       
00402     } else {
00403     
00404       // If we are less than mag 3 and negative, then we are suddenly
00405       // looking at decimal numbers not in scientific notation.
00406       // Therefore we hold as many decimal places as the magnitude.
00407       
00408       decimals = static_cast<int>( abs( rmag ) );
00409       
00410     }
00411   
00412   }
00413   // @todo decimals should never be negative here, but it does end up
00414   //    negative in some cases. See the "dirty fix" in Range.cxx, that
00415   //    dirty-fixed this problem too. But a better fix is needed. 
00416   if (decimals < 0) 
00417     decimals = 0;
00418    
00419   sprintf( pstr, "%%1.%df", decimals );
00420 
00421   y = first_tick;
00422   const Range & range = axis.getRange(false);
00423   double range_high = range.high();
00424   range_high *= scale_factor;
00425   
00426   while( y <= range_high || FLT_EQUAL( range_high, y ) )
00427     {
00428       double value = 0.0;
00429       
00430       if ( axistype == Axes::X ) {
00431         value = moduloAddX( y, xOffset() );
00432       }
00433       else if ( axistype == Axes::Y ) {
00434         value = moduloAddY( y, yOffset() );
00435       }
00436       
00437       // Now we either keep the original magnitude
00438       // or reduce it in order to express it in scientific notation.
00439       
00440       if ( use_pmag ) ylabel = value / pow( 10.0, pmag );
00441       else ylabel = value;
00442       
00443       value /= scale_factor;
00444       sprintf( labl, pstr, ylabel );
00445       m_ticks.push_back( AxisTick ( y, labl ) );
00446       
00447       num_ticks++;
00448       if ( tick_step == 0.0 ) break;
00449       y += tick_step;
00450       
00451     }
00452   
00453   return m_ticks;
00454 }
00455 
00456 
00457 bool
00458 PeriodicBinaryTransform::
00459 isLinearInXY () const
00460 {
00461   return false;
00462 }
00463 
00464 
00465 
00466 
00467 
00468 /* Use WCSLIB to do transform. */
00469 
00470 void 
00471 PeriodicBinaryTransform::
00472 initwcs(const std::string &transformName, double* crpix, 
00473         double* crval, double* cdelt, double crota2, bool galactic)
00474 {
00475 #ifdef HAVE_WCSLIB
00476     // Assert the sizeof wcsprm struct is smaller than the size allocated. 
00477     assert(sizeof(wcsprm)<=2000);
00478     m_wcs = reinterpret_cast<wcsprm*>(m_wcs_struct);
00479     m_wcs->flag = -1;
00480 
00481     // Call wcsini() in WCSLIB.
00482     int naxis = 2;
00483     wcsini(1, naxis, m_wcs);
00484 
00485     // Transfer parameters from c++ class to c stuct.
00486     std::string 
00487         lon_type = (galactic? "GLON-" : "RA---") + transformName,
00488         lat_type =  (galactic? "GLAT-" : "DEC--") + transformName;
00489     strcpy(m_wcs->ctype[0], lon_type.c_str() );
00490     strcpy(m_wcs->ctype[1], lat_type.c_str() );
00491 
00492     for( int i=0; i<naxis; ++i){
00493         m_wcs->crval[i] = crval[i]; // reference value
00494         m_wcs->crpix[i] = crpix[i]; // pixel coordinate
00495         m_wcs->cdelt[i] = cdelt[i]; // scale factor
00496     }
00497 #else
00498     throwWCSMissing ();
00499 #endif
00500 
00501 }
00502 
00503 void
00504 PeriodicBinaryTransform::
00505 throwWCSMissing () const
00506 {
00507   const std::string what
00508     ( "HippoDraw has not been compiled with WCSLIB support" );
00509   throw std::runtime_error ( what );
00510 }
00511 
00512 
00513 
00514 /* virtual */
00515 void
00516 PeriodicBinaryTransform::
00517 transform ( double & lon, double & lat ) const
00518 {
00519 
00520 #ifdef HAVE_WCSLIB
00521     int ncoords = 1;
00522     int nelem = 2;
00523     double  imgcrd[2], pixcrd[2];
00524     double phi[1], theta[1];
00525     int stat[1];
00526 
00527     double worldcrd[] ={lon,lat};
00528 
00529     //    int returncode = 
00530     wcss2p(m_wcs, ncoords, nelem, worldcrd, phi, theta, imgcrd, pixcrd, stat);
00531     
00532     lon = pixcrd[0];
00533     lat = pixcrd[1]; 
00534 #else
00535     throwWCSMissing ();
00536 #endif
00537 
00538 
00539 }
00540 
00541 
00542 bool
00543 PeriodicBinaryTransform::
00544 inverseTransform ( double & lon, double & lat ) const
00545 {
00546   
00547 #ifdef HAVE_WCSLIB
00548     int ncoords = 1;
00549     int nelem = 2;
00550     double worldcrd[2], imgcrd[2];
00551     double phi[1], theta[1];
00552     int stat[1];
00553 
00554     double pixcrd[] = {lon,lat};
00555 
00556     int returncode = 
00557     wcsp2s(m_wcs, ncoords, nelem, pixcrd, imgcrd, phi, theta, worldcrd, stat);
00558     
00559 
00560     // If there is invalid coordinate, return false.
00561     if ( returncode == 8 ) return false;
00562 
00563     // Make sure the inverse transform result in (-180,180) or (0,360)
00564     //if (worldcrd[0]<0) worldcrd[0]=worldcrd[0]+360;
00565     lon = worldcrd [0];
00566     lat = worldcrd [1];    
00567         
00568 
00569     return true;
00570 #else
00571     throwWCSMissing ();
00572     return false;
00573 #endif
00574 }
00575 
00576 void
00577 PeriodicBinaryTransform::
00578 transform ( std::vector< double > & lon,
00579             std::vector< double > & lat ) const
00580 {     
00581 
00582 #ifdef HAVE_WCSLIB
00583   assert ( lat.size() == lon.size() );
00584   int size = lat.size();
00585 
00586   // Form vector to array.
00587   vector < double >  worldcrd ( 2*size );
00588   for ( int i = 0; i < size; i++ ) {
00589     worldcrd[2*i]= lon[i];
00590     worldcrd[2*i+1]= lat[i];
00591   }
00592   
00593   int ncoords = size;
00594   int nelem = 2;
00595 
00596   vector < double > pixcrd( 2*size );
00597   vector < double > imgcrd( 2*size );
00598   vector < double > phi( size );
00599   vector < double > theta( size );
00600   vector < int > stat ( size );
00601 
00602   //int returncode =
00603   wcss2p ( m_wcs, ncoords, nelem, &worldcrd[0], &phi[0], &theta[0],
00604            &imgcrd[0], &pixcrd[0], &stat[0] );
00605   
00606   // From array to vector;  
00607   for ( int i = 0; i < size; i++ ) {
00608     lon[i]=pixcrd[2*i];
00609     lat[i]=pixcrd[2*i+1];
00610   }
00611    
00612 #else
00613     throwWCSMissing ();
00614 #endif
00615 }

Generated for HippoDraw Class Library by doxygen