PeriodicBinaryTransform.cxx
Go to the documentation of this file.
1 
11 // for wcslib
12 #ifdef HAVE_CONFIG_H
13 #include "config.h"
14 #endif
15 
17 
18 #include "axes/AxisModelBase.h"
19 #include "graphics/Rectangle.h"
20 #include "UnaryTransform.h"
21 
22 #ifdef HAVE_WCSLIB
23 #ifdef _WIN32
24 #include "wcslib/C/wcs.h"
25 #else
26 #include "wcslib/wcs.h"
27 #endif
28 #endif
29 
30 #include <cmath>
31 #include <cstdio>
32 #include <stdexcept>
33 #include <sstream>
34 #include <cassert>
35 
36 using namespace hippodraw;
37 
38 using std::max;
39 using std::abs;
40 using std::vector;
41 
49  m_x_limits ( -180.0, +180.0 ),
50  m_y_limits ( - 90.0, + 90.0 ),
51  m_x_offset ( 0.0 ),
52  m_y_offset ( 0.0 )
53 {
54 }
55 
57  bool is_periodic,
58  bool needs_grid,
59  bool needs_x_ticks,
60  bool needs_y_ticks,
61  double xlo, double xhi,
62  double ylo, double yhi) :
63  BinaryTransform ( z, is_periodic, needs_grid, needs_x_ticks, needs_y_ticks ),
64  m_x_limits ( xlo, xhi ),
65  m_y_limits ( ylo, yhi ),
66  m_x_offset ( 0.0 ),
67  m_y_offset ( 0.0 )
68 {
69 }
70 
73  BinaryTransform ( t ),
74  m_x_limits( t.limitX() ),
75  m_y_limits( t.limitY() ),
76  m_x_offset( t.xOffset() ),
77  m_y_offset( t.yOffset() )
78 {
79 }
80 
82 {
83  // Does nothing
84 }
85 
87 {
88  return m_x_limits;
89 }
90 
91 
93 {
94  return m_y_limits;
95 }
96 
97 
99 {
100  return m_x_offset;
101 }
102 
103 void PeriodicBinaryTransform::setXOffset( double x_offset )
104 {
105  m_x_offset = x_offset;
106 }
107 
109 {
110  return m_y_offset;
111 }
112 
113 void PeriodicBinaryTransform::setYOffset( double y_offset )
114 {
115 
116  m_y_offset = y_offset;
117 }
118 
119 
120 double
122 moduloAdd( double a1, double a2, hippodraw::Axes::Type axis ) const
123 {
124  if ( axis == Axes::X ) {
125  return moduloAddX ( a1, a2);
126  }
127  else if ( axis == Axes::Y ) {
128  return moduloAddY ( a1, a2);
129  }
130  else return a1 + a2;
131 }
132 
133 double
135 moduloSub( double s1, double s2, hippodraw::Axes::Type axis ) const
136 {
137  if ( axis == Axes::X ) {
138  return moduloSubX ( s1, s2);
139  }
140  else if ( axis == Axes::Y ) {
141  return moduloSubY ( s1, s2);
142  }
143  else return s1 - s2;
144 }
145 
146 double PeriodicBinaryTransform::moduloAddX( double x1, double x2 ) const
147 {
148  if ( x2 < -DBL_EPSILON )
149  return moduloSubX ( x1, -x2 );
150 
151  double overshoot = x1 + x2 - m_x_limits.high();
152 
153  if ( overshoot > DBL_EPSILON ) {
154  return m_x_limits.low() + overshoot;
155  }
156  else {
157  return x1 + x2;
158  }
159 }
160 
161 double PeriodicBinaryTransform::moduloSubX( double x1, double x2 ) const
162 {
163  if ( x2 < -DBL_EPSILON )
164  return moduloAddX ( x1, -x2 );
165 
166  double undershoot = m_x_limits.low() - ( x1 - x2 );
167 
168  if ( undershoot > DBL_EPSILON)
169  return m_x_limits.high() - undershoot;
170  else
171  return x1 - x2;
172 }
173 
174 
175 double PeriodicBinaryTransform::moduloAddY( double y1, double y2 ) const
176 {
177  if ( y2 < -DBL_EPSILON ) return moduloSubY ( y1, -y2 );
178 
179  double overshoot = y1 + y2 - m_y_limits.high();
180 
181  if ( overshoot > DBL_EPSILON ) {
182  return m_y_limits.low() + overshoot;
183  }
184  else {
185  return y1 + y2;
186  }
187 
188 }
189 
190 double PeriodicBinaryTransform::moduloSubY( double y1, double y2 ) const
191 {
192  if ( y2 < -DBL_EPSILON ) return moduloAddY ( y1, -y2 );
193 
194  double undershoot = m_y_limits.low() - ( y1 - y2 );
195 
196  if ( undershoot > DBL_EPSILON ) return m_y_limits.high() - undershoot;
197  else return y1 - y2;
198 }
199 
200 
202  const Range & lon )
203 {
204  double x_lo = lat.low ();
205  double x_hi = lat.high ();
206 
207  double y_lo = lon.low ();
208  double y_hi = lon.high ();
209 
210  double x, y;
211 
212  double x_min = 1000;
213  double x_max = -1000;
214  double y_min = 1000;
215  double y_max = -1000;
216 
217  int n = 50;
218  double dx = ( x_hi - x_lo ) / n;
219  double dy = ( y_hi - y_lo ) / n;
220 
221 
222  // Finding out xmin, xmax, ymin, ymax along line y = y_lo
223  for ( int i = 0; i <= n; i++)
224  {
225  x = x_lo + i * dx;
226  y = y_lo;
227  transform ( x, y );
228  x_min = ( x_min < x ) ? x_min : x;
229  x_max = ( x_max > x ) ? x_max : x;
230  y_min = ( y_min < y ) ? y_min : y;
231  y_max = ( y_max > y ) ? y_max : y;
232  }
233 
234  // Finding out xmin, xmax, ymin, ymax along line y = y_hi
235  for ( int i = 0; i <= n; i++)
236  {
237  x = x_lo + i * dx;
238  y = y_hi;
239  transform ( x, y );
240  x_min = ( x_min < x ) ? x_min : x;
241  x_max = ( x_max > x ) ? x_max : x;
242  y_min = ( y_min < y ) ? y_min : y;
243  y_max = ( y_max > y ) ? y_max : y;
244  }
245 
246  // Finding out xmin, xmax, ymin, ymax along line x = x_lo
247  for ( int i = 0; i <= n; i++)
248  {
249  x = x_lo;
250  y = y_lo + i * dy;
251  transform ( x, y );
252  x_min = ( x_min < x ) ? x_min : x;
253  x_max = ( x_max > x ) ? x_max : x;
254  y_min = ( y_min < y ) ? y_min : y;
255  y_max = ( y_max > y ) ? y_max : y;
256  }
257 
258  // Finding out xmin, xmax, ymin, ymax along line x = x_hi
259  for ( int i = 0; i <= n; i++)
260  {
261  x = x_hi;
262  y = y_lo + i * dy;
263  transform ( x, y );
264  x_min = ( x_min < x ) ? x_min : x;
265  x_max = ( x_max > x ) ? x_max : x;
266  y_min = ( y_min < y ) ? y_min : y;
267  y_max = ( y_max > y ) ? y_max : y;
268  }
269 
270  return Rect ( x_min, y_min, x_max - x_min, y_max - y_min );
271 }
272 
273 
274 void PeriodicBinaryTransform::validate ( Range & lat, Range & lon ) const
275 {
276  if ( lat.low () < m_x_limits.low() ) lat.setLow ( m_x_limits.low() );
277  if ( lat.high () > m_x_limits.high() ) lat.setHigh ( m_x_limits.high() );
278 
279  if ( lon.low () < m_y_limits.low() ) lon.setLow ( m_y_limits.low() );
280  if ( lon.high () > m_y_limits.high() ) lon.setHigh ( m_y_limits.high() );
281 }
282 
286 const vector < AxisTick > &
289 {
290  if ( axis == Axes::Z ) {
291  return m_z -> setTicks ( model );
292  }
293  //else
294  setTickStep ( model );
295  setFirstTick ( model );
296  return genTicks ( model, axis );
297 }
298 
300 void
304  const Range & )
305 {
306  // Does not do anything as of now
307  return;
308 }
309 
310 inline double FLT_EQUAL( double x, double y )
311 {
312  return ( (double)abs( x - y ) <= 2.0 * ( y * FLT_EPSILON + FLT_MIN ) );
313 }
314 
315 
317 {
318  const Range & range = axis.getRange(false);
319  double rangeLength = range.length();
320 
321  double scale_factor = axis.getScaleFactor();
322  rangeLength *= scale_factor;
323 
324  // The following algorithm determines the magnitude of the range...
325  double rmag = floor( log10( rangeLength ) );
326  axis.setRMag( rmag );
327 
328  double scalelow = range.low() * scale_factor;
329  double scalehigh = range.high() * scale_factor;
330 
331  // We will also need the largest magnitude for labels.
332  double pmag = max( floor( log10( abs ( scalehigh ) ) ),
333  floor( log10( abs ( scalelow ) ) ) );
334  axis.setPMag( pmag );
335 
336  axis.setTickStep( rangeLength / 4.0 );
337 }
338 
340 {
341  const Range & range = axis.getRange(false);
342 
343  axis.setFirstTick ( range.low() );
344 }
345 
346 
351 const vector < AxisTick > &
354 {
355  double y = 0.0, ylabel;
356 
357  int num_ticks = 0;
358  m_ticks.clear();
359  double pmag = axis.getPMag();
360  double rmag = axis.getRMag();
361  double first_tick = axis.getFirstTick();
362  double tick_step = axis.getTickStep();
363  double scale_factor = axis.getScaleFactor();
364 
365  // pmag will get set to 0 if it is less than or equal to 3. This
366  // is used later to determine scientific notation. However, m_rmag
367  // is still needed as the original magnitude for calculations such
368  // as decimal place notation, and rounding to nice numbers.
369 
370  bool use_pmag = abs ( pmag ) > 3.0;
371 
372  axis.setUsePMag ( use_pmag );
373 
374  char pstr[10];
375  char labl[10];
376 
377  int decimals = 0;
378 
379  // The following if-else block decides the pstr string, which holds
380  // the number of decimal places in the label.
381 
382  // if( fabs( m_pmag ) > 3.0 ) {
383  if ( use_pmag ) {
384  // If we are greater than mag 3, we are showing scientific
385  // notation. How many decimals we show is determined by the
386  // difference between the range magnitude and the power magnitude.
387 
388  decimals = static_cast<int>( pmag - rmag );
389  // minumum 1 decimal in scientific notation
390 
391  if( !decimals ) decimals++;
392 
393  } else {
394 
395  if( rmag > 0.0 ){
396 
397  // If we are less than mag 3 and positive, then no decimal
398  // accuracy is needed.
399 
400  decimals = 0;
401 
402  } else {
403 
404  // If we are less than mag 3 and negative, then we are suddenly
405  // looking at decimal numbers not in scientific notation.
406  // Therefore we hold as many decimal places as the magnitude.
407 
408  decimals = static_cast<int>( abs( rmag ) );
409 
410  }
411 
412  }
413  // @todo decimals should never be negative here, but it does end up
414  // negative in some cases. See the "dirty fix" in Range.cxx, that
415  // dirty-fixed this problem too. But a better fix is needed.
416  if (decimals < 0)
417  decimals = 0;
418 
419  sprintf( pstr, "%%1.%df", decimals );
420 
421  y = first_tick;
422  const Range & range = axis.getRange(false);
423  double range_high = range.high();
424  range_high *= scale_factor;
425 
426  while( y <= range_high || FLT_EQUAL( range_high, y ) )
427  {
428  double value = 0.0;
429 
430  if ( axistype == Axes::X ) {
431  value = moduloAddX( y, xOffset() );
432  }
433  else if ( axistype == Axes::Y ) {
434  value = moduloAddY( y, yOffset() );
435  }
436 
437  // Now we either keep the original magnitude
438  // or reduce it in order to express it in scientific notation.
439 
440  if ( use_pmag ) ylabel = value / pow( 10.0, pmag );
441  else ylabel = value;
442 
443  value /= scale_factor;
444  sprintf( labl, pstr, ylabel );
445  m_ticks.push_back( AxisTick ( y, labl ) );
446 
447  num_ticks++;
448  if ( tick_step == 0.0 ) break;
449  y += tick_step;
450 
451  }
452 
453  return m_ticks;
454 }
455 
456 
457 bool
459 isLinearInXY () const
460 {
461  return false;
462 }
463 
464 
465 
466 
467 
468 /* Use WCSLIB to do transform. */
469 
470 void
472 initwcs(const std::string &transformName, double* crpix,
473  double* crval, double* cdelt, double crota2, bool galactic)
474 {
475 #ifdef HAVE_WCSLIB
476  // Assert the sizeof wcsprm struct is smaller than the size allocated.
477  assert(sizeof(wcsprm)<=2000);
478  m_wcs = reinterpret_cast<wcsprm*>(m_wcs_struct);
479  m_wcs->flag = -1;
480 
481  // Call wcsini() in WCSLIB.
482  int naxis = 2;
483  wcsini(1, naxis, m_wcs);
484 
485  // Transfer parameters from c++ class to c stuct.
486  std::string
487  lon_type = (galactic? "GLON-" : "RA---") + transformName,
488  lat_type = (galactic? "GLAT-" : "DEC--") + transformName;
489  strcpy(m_wcs->ctype[0], lon_type.c_str() );
490  strcpy(m_wcs->ctype[1], lat_type.c_str() );
491 
492  for( int i=0; i<naxis; ++i){
493  m_wcs->crval[i] = crval[i]; // reference value
494  m_wcs->crpix[i] = crpix[i]; // pixel coordinate
495  m_wcs->cdelt[i] = cdelt[i]; // scale factor
496  }
497 #else
498  throwWCSMissing ();
499 #endif
500 
501 }
502 
503 void
506 {
507  const std::string what
508  ( "HippoDraw has not been compiled with WCSLIB support" );
509  throw std::runtime_error ( what );
510 }
511 
512 
513 
514 /* virtual */
515 void
517 transform ( double & lon, double & lat ) const
518 {
519 
520 #ifdef HAVE_WCSLIB
521  int ncoords = 1;
522  int nelem = 2;
523  double imgcrd[2], pixcrd[2];
524  double phi[1], theta[1];
525  int stat[1];
526 
527  double worldcrd[] ={lon,lat};
528 
529  // int returncode =
530  wcss2p(m_wcs, ncoords, nelem, worldcrd, phi, theta, imgcrd, pixcrd, stat);
531 
532  lon = pixcrd[0];
533  lat = pixcrd[1];
534 #else
535  throwWCSMissing ();
536 #endif
537 
538 
539 }
540 
541 
542 bool
544 inverseTransform ( double & lon, double & lat ) const
545 {
546 
547 #ifdef HAVE_WCSLIB
548  int ncoords = 1;
549  int nelem = 2;
550  double worldcrd[2], imgcrd[2];
551  double phi[1], theta[1];
552  int stat[1];
553 
554  double pixcrd[] = {lon,lat};
555 
556  int returncode =
557  wcsp2s(m_wcs, ncoords, nelem, pixcrd, imgcrd, phi, theta, worldcrd, stat);
558 
559 
560  // If there is invalid coordinate, return false.
561  if ( returncode == 8 ) return false;
562 
563  // Make sure the inverse transform result in (-180,180) or (0,360)
564  //if (worldcrd[0]<0) worldcrd[0]=worldcrd[0]+360;
565  lon = worldcrd [0];
566  lat = worldcrd [1];
567 
568 
569  return true;
570 #else
571  throwWCSMissing ();
572  return false;
573 #endif
574 }
575 
576 void
578 transform ( std::vector< double > & lon,
579  std::vector< double > & lat ) const
580 {
581 
582 #ifdef HAVE_WCSLIB
583  assert ( lat.size() == lon.size() );
584  int size = lat.size();
585 
586  // Form vector to array.
587  vector < double > worldcrd ( 2*size );
588  for ( int i = 0; i < size; i++ ) {
589  worldcrd[2*i]= lon[i];
590  worldcrd[2*i+1]= lat[i];
591  }
592 
593  int ncoords = size;
594  int nelem = 2;
595 
596  vector < double > pixcrd( 2*size );
597  vector < double > imgcrd( 2*size );
598  vector < double > phi( size );
599  vector < double > theta( size );
600  vector < int > stat ( size );
601 
602  //int returncode =
603  wcss2p ( m_wcs, ncoords, nelem, &worldcrd[0], &phi[0], &theta[0],
604  &imgcrd[0], &pixcrd[0], &stat[0] );
605 
606  // From array to vector;
607  for ( int i = 0; i < size; i++ ) {
608  lon[i]=pixcrd[2*i];
609  lat[i]=pixcrd[2*i+1];
610  }
611 
612 #else
613  throwWCSMissing ();
614 #endif
615 }

Generated for HippoDraw Class Library by doxygen