FitsFile.cxx
Go to the documentation of this file.
1 
12 #include "FitsFile.h"
13 
14 #include "pattern/string_convert.h"
15 
16 #include <cstring>
17 #include <algorithm>
18 #include <stdexcept>
19 
20 #include <cassert>
21 
22 using std::string;
23 using std::vector;
24 
25 using namespace hippodraw;
26 
27 FitsFile::FitsFile ( const std::string & filename, bool write )
28  : FitsFileBase ( filename, write )
29 {
30 }
31 
32 void
34 fillHDUNames ( std::vector < std::string > & names )
35 {
36  names.clear ();
37  int number = getNumberOfHDU ( );
38 
39  for ( int i = 0; i < number; i++ ) {
40  moveToHDU ( i );
41  string value = stringValueForKey ( "EXTNAME" );
42 
43  if ( value == "" ) {
44  if ( i == 0 ) {
45  int number = intValueForKey ( "NAXIS" );
46  if ( number == 0 ) {
47  value = "Empty image";
48  }
49  }
50  }
51 
52  if ( value == "" ) value = "<no name>";
53  names.push_back ( value );
54  }
55 }
56 
57 int
59 fillColumnNames ( std::vector < std::string > & labels )
60 {
61  HduType type = getHduType ();
62  if ( type == Image ) {
64  }
65  else {
67  }
68 
69  return m_status;
70 }
71 
72 int
74 fillColumnNamesFromTable ( std::vector < std::string > & labels )
75 {
76  labels.clear ();
77  int columns = getNumberOfColumns ();
78 
79  if ( columns != 0 ) {
80  m_status = 0; // must NOT be reset on each fits call, else matchs
81  // same column on each call
82  for ( int i = 0; i < columns; i++ ) { // count like Fortran
83  char colname [ FLEN_VALUE ];
84  int col;
85  fits_get_colname ( m_fptr, CASEINSEN, const_cast< char * > ("*"),
86  colname, &col, &m_status );
87  labels.push_back ( string ( colname ) );
88  }
89  }
90  else { // must be image
91  labels.push_back ( string ( "none" ) );
92  }
93 
94  return m_status;
95 }
96 
97 int
99 fillColumnNamesFromImage ( std::vector < std::string > & labels )
100 {
101  labels.clear ();
102  int dimension = getImageDimensions ();
103  if ( dimension == 2 ) {
104  labels.push_back ( "Values" );
105  }
106  else if ( dimension == 3 ) {
107  const char type[] = "CTYPE3";
108  bool yes = hasKey ( type );
109  if ( yes ) {
110  const string name = stringValueForKey ( type );
111  int number = intValueForKey ( "NAXIS3" );
112  for ( int i = 0; i < number; i++ ) {
113  const string iv = String::convert ( i );
114  labels.push_back ( name + " " + iv );
115  }
116  }
117  }
118 
119  return m_status;
120 }
121 
122 int
124 fillDoubleVectorFromColumn ( std::vector < double > & vec,
125  int column )
126 {
127  int retval = 0;
128 
129  HduType type = getHduType ();
130  if ( type == Btable || type == Atable ) {
131  retval = fillFromTableColumn ( vec, column );
132  }
133  else {
134  retval = fillFromImage ( vec, column );
135  }
136 
137  return retval;
138 }
139 
140 int
142 fillFromTableColumn ( std::vector < double > & vec, int column )
143 {
144  int anynul;
145  double nulval = 0;
146  long nelements = vec.size();
147  vector<double>::iterator it = vec.begin();
148  double * ptr = &*it;
149  // Careful, counting columns like Fortran
150  m_status = 0;
151  fits_read_col_dbl( m_fptr, column+1, 1, 1, nelements, nulval,
152  ptr, &anynul, &m_status);
153 
154  return m_status;
155 }
156 
157 void
159 fillShape ( std::vector < intptr_t > & shape, int column )
160 {
161  m_status = 0; // must always set
162 
163  string tdimn ( "TDIM" );
164  tdimn += String::convert ( column+1 ); // a la Fortran
165  char value[80]; // max length of keyword record
166  fits_read_keyword ( m_fptr,
167  const_cast < char * > ( tdimn.c_str() ),
168  value,
169  NULL, & m_status );
170 
171  long rows = getNumberOfRows ();
172  if ( value[0] != 0 ) {
173  const string dims = value;
174  unsigned int count = std::count ( dims.begin(), dims.end(), ',' );
175  count += 1; // one for simple vector
176  int icount = 0;
177  vector < long > naxes ( count );
178  fits_read_tdim ( m_fptr, column+1, count, // input
179  &icount, &naxes[0], & m_status ); //output
180  shape.resize ( count + 1 );
181  shape [0] = rows;
182  std::copy ( naxes.begin(), naxes.end(), shape.begin() + 1 );
183  }
184  else { // no array
185  shape.resize ( 1, rows );
186  }
187 }
188 
189 int
191 fillAxisSizes ( std::vector < long > & vec ) const
192 {
193  vec.clear ();
194  int naxis = getImageDimensions ();
195  vec.resize ( naxis );
196 
197  vector < long > ::iterator first = vec.begin();
198  long * ptr = & *first;
199  m_status = 0;
200  fits_get_img_size ( m_fptr, naxis, ptr, & m_status );
201  assert ( m_status == 0 );
202 
203  return m_status;
204 }
205 
206 void
208 fillImageDeltas ( std::vector < double > & deltas ) const
209 {
210  deltas.clear ();
211  int naxis = getImageDimensions ();
212  deltas.reserve ( naxis );
213 
214  char key [ FLEN_KEYWORD];
215  char * keyroot = const_cast < char * > ("CDELT");
216  for ( int i = 0; i < naxis; i++ ) {
217  m_status = 0;
218  fits_make_keyn ( keyroot, i+1, key, & m_status );
219  bool yes = hasKey ( key );
220 
221  if ( yes ) {
222  double value = doubleValueForKey ( key );
223  deltas.push_back ( value );
224  }
225  else { // no key, take default
226  deltas.push_back ( 1.0 );
227  }
228  }
229 }
230 
231 void
233 fillRefPixelIndices ( std::vector < int > & indices ) const
234 {
235  indices.clear ();
236  int naxis = getImageDimensions ();
237  indices.reserve ( naxis );
238 
239  char key [ FLEN_KEYWORD];
240  char * keyroot = const_cast < char * > ( "CRPIX" );
241  for ( int i = 0; i < naxis; i++ ) {
242  m_status = 0;
243  fits_make_keyn ( keyroot, i+1, key, & m_status );
244  bool yes = hasKey ( key );
245 
246  if ( yes ) {
247  int value = intValueForKey ( key );
248  indices.push_back ( value );
249  }
250  else { // no key, take default
251  indices.push_back ( 1 );
252  }
253  }
254 }
255 
256 void
258 fillRefPixelValues ( std::vector < double > & values ) const
259 {
260  values.clear ();
261  int naxis = getImageDimensions ();
262  values.reserve ( naxis );
263 
264  char key [ FLEN_KEYWORD];
265  char * keyroot = const_cast < char * > ( "CRVAL" );
266  for ( int i = 0; i < naxis; i++ ) {
267  m_status = 0;
268  fits_make_keyn ( keyroot, i+1, key, & m_status );
269  bool yes = hasKey ( key );
270 
271  if ( yes ) {
272  double value = doubleValueForKey ( key );
273  values.push_back ( value );
274  }
275  else { // no key, take default
276  values.push_back ( 0. );
277  }
278  }
279 }
280 
281 bool
284 {
285  bool hammer = true;
286 
287  char key [FLEN_KEYWORD];
288  char * keyroot = const_cast < char * > ( "CTYPE" );
289  for ( int i = 0; i < 2; i++ ) {
290  fits_make_keyn ( keyroot, i+1, key, & m_status );
291  bool yes = hasKey ( key );
292 
293  if ( yes ) {
294  string value = stringValueForKey ( key );
295  if ( value.find ( "-AIT" ) == string::npos ) {
296  hammer = false;
297  }
298  } else {
299  hammer = false;
300  }
301  }
302 
303  return hammer;
304 }
305 
306 int
308 fillFromImage ( std::vector < double > & vec, unsigned int zplane )
309 {
310  vector < long > naxes;
311  fillAxisSizes ( naxes );
312 
313  int datatype = Double;
314  long nelements = naxes[0] * naxes[1];
315  double nulval = 0;
316  vector < double >::iterator first = vec.begin();
317  double * ptr = & *first;
318  int anynul;
319  m_status = 0;
320  if ( naxes.size () == 2 ) {
321  long fpixel[2] = { 1, 1 }; // index like Fortran
322  fits_read_pix ( m_fptr, datatype, fpixel, nelements,
323  & nulval, ptr, & anynul, & m_status );
324  } else {
325  long fpixel[] = { 1, 1, 1 }; // index like Fortran
326  fpixel[2] = zplane + 1; // like Fortran
327  fits_read_pix ( m_fptr, datatype, fpixel, nelements,
328  & nulval, ptr, & anynul, & m_status );
329  }
330 
331  return m_status;
332 }
333 
334 int
336 fillIntVectorFromColumn( std::vector < int > & vec, int column )
337 {
338  int anynul, status = 0;
339  int nulval = 0;
340  long nelements = vec.size();
341  vector<int>::iterator it = vec.begin();
342  int * ptr = new int [ nelements ];
343  // Careful, counting columns like Fortran
344  m_status = 0;
345  fits_read_col_int( m_fptr, column+1, 1, 1, nelements, nulval,
346  ptr, &anynul, &status);
347  copy ( ptr, ptr+nelements, it );
348 
349  return status;
350 }
351 
352 void
354 writeHDU ( long rows, int columns,
355  const std::vector < std::string > & names,
356  const std::vector < std::vector < int > > & shapes,
357  const std::string & extname )
358 {
359  char ** types = new char * [ columns ];
360  char ** tform = new char * [ columns ];
361  char ** tunits = 0;
362  char * m_extname;
363 
364  vector < string > forms;
365  for ( int i = 0; i < columns; i ++ ) {
366  types[i]=const_cast<char*> (names[i].c_str());
367  int size = 1;
368  const vector < int > & shape = shapes [ i ];
369  unsigned int rank = shape.size();
370  if ( rank > 1 ) {
371  for ( unsigned int j = 1; j < rank; j++ ) { // strip row index
372  int dim = shape [ j ];
373  size *= dim;
374  }
375  }
376  string form = String::convert ( size );
377  form += "D";
378  forms.push_back ( form );
379 // unsigned int n = form.size();
380  tform[i] = strdup ( form.c_str() );
381 
382 // tform [i] = const_cast <char *> ( forms[i].c_str() ); // will stay in scope
383  }
384 
385  m_extname = const_cast<char*> (extname.c_str());
386 
387  fits_create_tbl( m_fptr, BINARY_TBL, 0, columns,
388  types, tform, tunits, m_extname, &m_status);
389 //fits_report_error ( stdout, m_status );
390  for ( int i = 0; i < columns; i ++ ) {
391  const vector < int > & shape = shapes [ i ];
392  unsigned int rank = shape.size();
393  if ( rank > 1 ) {
394  vector < long > naxes;
395  for ( unsigned int j = 1; j < rank; j++ ) { // strip row index
396  naxes.push_back ( shape [ j ] );
397  }
398  fits_write_tdim ( m_fptr, i+1, // like Fortran
399  rank -1, &naxes [0], &m_status );
400  }
401  }
402  // clean up
403  delete types;
404  delete tform;
405 }
406 
407 void
409 writeImageHDU ( long x, long y )
410 {
411  long naxes[2]={x,y};
412  fits_create_img ( m_fptr, DOUBLE_IMG, 2, naxes, &m_status );
413 }
414 
415 
416 void
418 writeColumn ( int col, const std::vector < double > & data )
419 {
420  LONGLONG nelements = data.size ();
421 LONGLONG firstrow = 1;
422 LONGLONG firstelem = 1;
423  fits_write_col( m_fptr, TDOUBLE, col+1, firstrow, firstelem, nelements,
424  const_cast < double * > ( &data[0] ), &m_status);
425 }
426 
427 void
429 writePix ( long x, long y, const std::vector < double > & data )
430 {
431  long fpixel[2]={1,1};
432  long nelements = x * y;
433  fits_write_pix ( m_fptr, TDOUBLE, fpixel, nelements,
434  const_cast < double * > ( &data[0] ), &m_status );
435 }
436 
437 void
440 {
441  fits_close_file( m_fptr, &m_status );
442  m_fptr = 0; // to remember its already closed.
443 }
444 
445 bool
447 pixCenter () const
448 {
449 
450  char * key = const_cast < char * > ( "PIXCENT" );
451  bool yes = hasKey ( key );
452  if ( yes )
453  {
454  if ( stringValueForKey ( key ) == "T" ) return true;
455  }
456 
457  return false;
458 
459 }
460 
461 void
463 writeRefPixelValues ( double value1, double value2 )
464 {
465  char * key1 = const_cast < char * > ( "CRVAL1" );
466  char * key2 = const_cast < char * > ( "CRVAL2" );
467 
468  fits_write_key( m_fptr, TDOUBLE, key1, &value1, NULL, &m_status );
469  fits_write_key( m_fptr, TDOUBLE, key2, &value2, NULL, &m_status );
470 }

Generated for HippoDraw Class Library by doxygen