FitsController.cxx
Go to the documentation of this file.
1 
12 #include "FitsController.h"
13 
14 #include "FitsFile.h"
15 #include "FitsNTuple.h"
16 
19 #include "datasrcs/TupleCut.h"
20 #include "pattern/string_convert.h"
21 #include "plotters/PlotterBase.h"
22 
23 #include <algorithm>
24 #include <fstream>
25 
26 
27 #include <cassert>
28 
29 using std::ifstream;
30 using std::runtime_error;
31 using std::string;
32 using std::vector;
33 
34 using namespace hippodraw;
35 
38 typedef std::map < std::string, FitsFile * >::iterator FileMapIterator;
39 
41 
45 {
46  if ( s_instance == 0 ) {
47  s_instance = new FitsController ();
48  }
49  return s_instance;
50 }
51 
54 {
55 }
56 
57 const std::string &
59 version () const
60 {
61  float version = 0.0;
62  fits_get_version ( & version );
63 
64  m_version = String::convert ( version );
65 
66  return m_version;
67 }
68 
69 FitsFile *
71 openFile ( const std::string & name )
72 {
73  FitsFile * file = 0;
74  FileMapIterator first = m_file_map.find ( name );
75 
76  if ( first == m_file_map.end() ) {
77  ifstream test ( name.c_str (), std::ios::in );
78  if ( test.is_open () == false ) {
79  string what ( "FitsController: File `" );
80  what += name;
81  what += "' not found";
82  throw std::runtime_error ( what );
83  }
84 
85  file = new FitsFile ( name );
86  int status = file -> status();
87 
88  if ( status != 0 ) {
89  string what ( "cfitsio: Error opening file `" );
90  what += name;
91  what += "'";
92  delete file;
93  throw runtime_error ( what );
94  }
95  m_file_map [ name ] = file;
96  }
97  else {
98  file = first -> second;
99  }
100 
101  return file;
102 }
103 
106 void
108 closeFile ( const std::string & )
109 {
110 // FileMapIterator where = m_file_map.find ( name );
111 // if ( where != m_file_map.end () ) {
112 // m_file_map.erase ( where );
113 // fitsfile * file = where -> second;
114 // fits_close_file ( file, & status );
115 // }
116 }
117 
118 const vector < string > &
120 getNTupleNames ( const std::string & file_name )
121 {
122  m_ntuple_names.clear ();
123  FitsFile * file = openFile ( file_name );
124  if ( file != 0 ) {
125  file -> fillHDUNames ( m_ntuple_names );
126  }
127 
128  return m_ntuple_names;
129 }
130 
131 
132 DataSource *
134 createNTuple ( const std::string & filename,
135  const std::string & name )
136 {
137  int index = -1;
138 
139  const vector < string > & names
140  = getNTupleNames ( filename );
141  unsigned int size = names.size ();
142 
143  for ( unsigned int i = 0; i < size; i++ ) {
144  if ( name == names[i] ) {
145  index = i;
146  break;
147  }
148  }
149 
150  if ( index == -1 ) {
151  string what ( "FitsController: File `" );
152  what += filename;
153  what += "' ";
154  what += "does not have extensison with name `";
155  what += name;
156  what += "'";
157  throw runtime_error ( what );
158  }
159 
160  return createNTuple ( filename, name, index );
161 }
162 
163 DataSource *
165 createNTuple ( const std::string & filename,
166  const std::string & name,
167  int index )
168 {
169  DataSource * ntuple = 0;
170 
171  string ds_name = filename;
172  ds_name += ": ";
173  ds_name += name;
174 
175  FitsFile * file = openFile ( filename );
176  if ( file == NULL ) return NULL;
177 
178  int retval = file -> moveToHDU ( index );
179  if ( retval == 0 ) {
180  try {
181  ntuple = new FitsNTuple ( file );
182  }
183  catch ( const runtime_error & e ) {
184  throw e;
185  }
186 
187  ntuple -> setTitle ( name );
188  ntuple -> setName ( ds_name );
190  controller -> registerNTuple ( ds_name, ntuple );
191  controller -> registerDataSourceFile ( ntuple );
192  }
193 
194  return ntuple;
195 }
196 
197 void
199 checkForImage ( PlotterBase * plotter, const DataSource & source )
200 {
201  const FitsFile * file = 0;
202  try {
203  const FitsNTuple & ntuple
204  = dynamic_cast < const FitsNTuple & > ( source );
205  file = ntuple.getFile ();
206  }
207  catch ( ... ) {
208  // do nothing
209  }
210  if ( file == 0 ) return;
211 
212  FitsFileBase::HduType type = file -> getHduType ();
213  if ( type != FitsFileBase::Image ) return;
214 
215  vector < long > sizes;
216  file -> fillAxisSizes ( sizes );
217  if ( sizes.size () > 3 ) {
218  string what ( "FitsController: greater then 3D images not supported" );
219  throw std::runtime_error ( what );
220  }
221 
222  plotter -> setNumberOfBins ( Axes::X, sizes[0] );
223  plotter -> setNumberOfBins ( Axes::Y, sizes[1] );
224 
225  vector < double > deltas;
226  file -> fillImageDeltas ( deltas );
227  plotter -> setBinWidth ( Axes::X, deltas[0] );
228  plotter -> setBinWidth ( Axes::Y, deltas[1] );
229 
230  vector < double > ref_values;
231  file -> fillRefPixelValues ( ref_values );
232  vector < int > ref_indices;
233  file -> fillRefPixelIndices ( ref_indices );
234 
235  // Process fits file with PIXCENT keyword.
236  double pixToRef;
237  if ((file ->pixCenter ())==true) pixToRef=0.5;
238  else pixToRef=0.0;
239 
240  double x_orig = - deltas[0] * ( ref_indices[0]-1+pixToRef ) + ref_values[0];
241  double y_orig = - deltas[1] * ( ref_indices[1]-1+pixToRef ) + ref_values[1];
242  plotter ->setOffset ( Axes::X, x_orig );
243  plotter ->setOffset ( Axes::Y, y_orig );
244 
245 
246  // Fix out of range bug.
247  Range range;
248  if ( deltas[0] < 0. ) {
249  range.setRange ( x_orig + sizes[0] * deltas[0], x_orig, -deltas[0] );
250  }
251  else {
252  range.setRange ( x_orig, x_orig + (sizes[0] ) * deltas[0], 1. );
253  }
254  plotter -> setRange ( Axes::X, range, false, true );
255  range.setLow ( y_orig );
256  range.setLength ( ( sizes[1] ) * deltas[1] );
257  plotter -> setRange ( Axes::Y, range, false, true );
258 
259  plotter -> matrixTranspose ( true );
260 
261  bool yes = file -> isHammerAitoff ();
262  if ( yes ) {
263  plotter ->setAspectRatio ( 2.0 );
264  plotter ->setFitsTransform ("HammerAito");
265  }
266 }
267 
270 void
272 writeNTupleToFile ( const DataSource * ntuple,
273  const std::string & filename )
274 {
275  FitsFile * file = new FitsFile ( filename, true );
276 
277  long row = static_cast<long> ( ntuple -> rows() );
278  int column = static_cast<int> (ntuple -> columns());
279  const vector <string> labels = ntuple -> getLabels();
280  vector < vector < int > > shapes;
281 
282  for ( int i = 0; i < column; i++ ){
283  vector < int > shape;
284  ntuple -> fillShape ( shape, i );
285  shapes.push_back ( shape );
286  }
287  string extname;
288  const string fullname = ntuple -> getName();
289  string::size_type pos1 = fullname.find_last_of ( ":" );
290  string::size_type pos2 = fullname.find_last_of ( "/" );
291 
292  if ( pos1 != string::npos && pos1!=1 ) {
293  extname = fullname.substr ( pos1+2 );
294  }
295  else if (pos2 != string::npos) { // only filename
296  extname = "<no name>";
297  }
298  else { // only title
299  extname = fullname;
300  }
301 
302  file -> writeHDU ( row, column, labels, shapes, extname );
303  for ( int i = 0; i < column; i++ ) {
304  const vector < double > & col = ntuple -> getColumn (i );
305  file -> writeColumn ( i, col );
306  }
307 
308  file -> writeCloseFile ();
309 
310  delete file;
311 
312 }
313 
314 void
316 writeNTupleToFile ( const std::string & name, const std::string & filename )
317 {
319  DataSource * ntuple
320  = controller -> findDataSource ( name ); // throws exception if not found
321  if ( ntuple == 0 ) return;
322 
323  writeNTupleToFile ( ntuple, filename );
324 }
325 
326 void
328 writeImageToFile ( unsigned int x, unsigned int y,
329  const std::vector <double> & data,
330  const std::string & filename )
331 {
332 
333  FitsFile * file = new FitsFile ( filename, true );
334 
335  long xx = static_cast<long> ( x );
336  long yy = static_cast<long> ( y );
337 
338  double xlow = 0.0 - static_cast<double> (x) / 2 ;
339  double ylow = 0.0 - static_cast<double> (y) / 2;
340 
341  file -> writeImageHDU ( xx,yy );
342  file -> writeRefPixelValues ( xlow, ylow );
343  file -> writePix ( xx, yy, data );
344 
345  file -> writeCloseFile ();
346 
347  delete file;
348 
349 }
350 
351 size_t
353 calcColumnWidth ( const DataSource * source, unsigned int column ) const
354 {
355  vector < int > shape;
356  source -> fillShape ( shape, column );
357  size_t rank = shape.size ();
358  size_t count = 1;
359 
360  for ( size_t i = 1; i < rank; i++ ) {
361  count *= shape[i];
362  }
363  return count;
364 }
365 
366 int
369  const std::string & filename,
370  const std::string & dsname,
371  const std::vector < std::string > & column_list,
372  const std::vector < const TupleCut * > & cut_list)
373 {
374  if ( column_list.empty() ) return 1;
375 
376  FitsFile * file = new FitsFile ( filename, true );
377 
378  unsigned int columnNumber = column_list.size();
379  unsigned int size = ds->rows();
380 
381  vector < bool > acceptArray;
382  CutController::fillAcceptedRows ( acceptArray, ds, cut_list );
383 
384  long rows = std::count ( acceptArray.begin(), acceptArray.end(), true );
385  int columns = static_cast<int> (columnNumber);
386 
387  // Following collected to create FITS table before writing.
388 
389  vector < int > col_indices ( columnNumber );
390 
391  for ( unsigned int i = 0; i < columnNumber; i++ ) {
392  const string & label = column_list [ i ];
393  int index = ds -> indexOf ( label );
394  if ( index < 0 ) {
395  ds -> throwIfInvalidLabel ( label );
396  }
397  col_indices [i] = index;
398  }
399 
400  vector < size_t > sizes;
401  vector < vector < int > > shapes;
402  for ( int i = 0; i < columns; i++ ) {
403  unsigned int size = calcColumnWidth ( ds, col_indices[i] );
404  sizes.push_back ( size );
405  vector < int > shape;
406  ds -> fillShape ( shape, col_indices[i] );
407  shapes.push_back ( shape );
408  }
409 
410  file -> writeHDU ( rows, columns, column_list,
411  shapes, dsname );
412 
413 
414  for ( unsigned int k = 0; k < columnNumber; k++ ) {
415  std::vector < double > tempColumn;
416  if ( sizes [ k ] == 1 ) {
417  tempColumn.reserve (size );
418  }
419  else {
420  tempColumn.reserve ( size * sizes[k] );
421  }
422  for ( unsigned int i = 0; i<size; i++ ) {
423  if ( acceptArray[i]==true ) {
424  int index = col_indices[k];
425  if ( sizes [ k ] == 1 ) {
426  tempColumn.push_back ( ds -> valueAtNoCache (i, index) );
427  }
428  else { // array variable
429  double * array = ds -> doubleArrayAt ( i, index );
430  for ( std::size_t l = 0; l < sizes [ k ]; l++ ) {
431  tempColumn.push_back ( array[l] );
432  }
433  }
434  }
435  }
436  file -> writeColumn ( k, tempColumn );
437 
438  tempColumn.clear();
439  }
440 
441  file -> writeCloseFile ();
442 
443  delete file;
444 
445  return 0;
446 }

Generated for HippoDraw Class Library by doxygen