00001
00011
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
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
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
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
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
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
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
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
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
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
00366
00367
00368
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
00380
00381
00382
00383 if ( use_pmag ) {
00384
00385
00386
00387
00388 decimals = static_cast<int>( pmag - rmag );
00389
00390
00391 if( !decimals ) decimals++;
00392
00393 } else {
00394
00395 if( rmag > 0.0 ){
00396
00397
00398
00399
00400 decimals = 0;
00401
00402 } else {
00403
00404
00405
00406
00407
00408 decimals = static_cast<int>( abs( rmag ) );
00409
00410 }
00411
00412 }
00413
00414
00415
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
00438
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
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
00477 assert(sizeof(wcsprm)<=2000);
00478 m_wcs = reinterpret_cast<wcsprm*>(m_wcs_struct);
00479 m_wcs->flag = -1;
00480
00481
00482 int naxis = 2;
00483 wcsini(1, naxis, m_wcs);
00484
00485
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];
00494 m_wcs->crpix[i] = crpix[i];
00495 m_wcs->cdelt[i] = cdelt[i];
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
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
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
00561 if ( returncode == 8 ) return false;
00562
00563
00564
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
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
00603 wcss2p ( m_wcs, ncoords, nelem, &worldcrd[0], &phi[0], &theta[0],
00604 &imgcrd[0], &pixcrd[0], &stat[0] );
00605
00606
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 }