00001
00012 #include "FitsFile.h"
00013
00014 #include "pattern/string_convert.h"
00015
00016 #include <cstring>
00017 #include <algorithm>
00018 #include <stdexcept>
00019
00020 #include <cassert>
00021
00022 using std::string;
00023 using std::vector;
00024
00025 using namespace hippodraw;
00026
00027 FitsFile::FitsFile ( const std::string & filename, bool write )
00028 : FitsFileBase ( filename, write )
00029 {
00030 }
00031
00032 void
00033 FitsFile::
00034 fillHDUNames ( std::vector < std::string > & names )
00035 {
00036 names.clear ();
00037 int number = getNumberOfHDU ( );
00038
00039 for ( int i = 0; i < number; i++ ) {
00040 moveToHDU ( i );
00041 string value = stringValueForKey ( "EXTNAME" );
00042
00043 if ( value == "" ) {
00044 if ( i == 0 ) {
00045 int number = intValueForKey ( "NAXIS" );
00046 if ( number == 0 ) {
00047 value = "Empty image";
00048 }
00049 }
00050 }
00051
00052 if ( value == "" ) value = "<no name>";
00053 names.push_back ( value );
00054 }
00055 }
00056
00057 int
00058 FitsFile::
00059 fillColumnNames ( std::vector < std::string > & labels )
00060 {
00061 HduType type = getHduType ();
00062 if ( type == Image ) {
00063 m_status = fillColumnNamesFromImage ( labels );
00064 }
00065 else {
00066 m_status = fillColumnNamesFromTable ( labels );
00067 }
00068
00069 return m_status;
00070 }
00071
00072 int
00073 FitsFile::
00074 fillColumnNamesFromTable ( std::vector < std::string > & labels )
00075 {
00076 labels.clear ();
00077 int columns = getNumberOfColumns ();
00078
00079 if ( columns != 0 ) {
00080 m_status = 0;
00081
00082 for ( int i = 0; i < columns; i++ ) {
00083 char colname [ FLEN_VALUE ];
00084 int col;
00085 fits_get_colname ( m_fptr, CASEINSEN, const_cast< char * > ("*"),
00086 colname, &col, &m_status );
00087 labels.push_back ( string ( colname ) );
00088 }
00089 }
00090 else {
00091 labels.push_back ( string ( "none" ) );
00092 }
00093
00094 return m_status;
00095 }
00096
00097 int
00098 FitsFile::
00099 fillColumnNamesFromImage ( std::vector < std::string > & labels )
00100 {
00101 labels.clear ();
00102 int dimension = getImageDimensions ();
00103 if ( dimension == 2 ) {
00104 labels.push_back ( "Values" );
00105 }
00106 else if ( dimension == 3 ) {
00107 const char type[] = "CTYPE3";
00108 bool yes = hasKey ( type );
00109 if ( yes ) {
00110 const string name = stringValueForKey ( type );
00111 int number = intValueForKey ( "NAXIS3" );
00112 for ( int i = 0; i < number; i++ ) {
00113 const string iv = String::convert ( i );
00114 labels.push_back ( name + " " + iv );
00115 }
00116 }
00117 }
00118
00119 return m_status;
00120 }
00121
00122 int
00123 FitsFile::
00124 fillDoubleVectorFromColumn ( std::vector < double > & vec,
00125 int column )
00126 {
00127 int retval = 0;
00128
00129 HduType type = getHduType ();
00130 if ( type == Btable || type == Atable ) {
00131 retval = fillFromTableColumn ( vec, column );
00132 }
00133 else {
00134 retval = fillFromImage ( vec, column );
00135 }
00136
00137 return retval;
00138 }
00139
00140 int
00141 FitsFile::
00142 fillFromTableColumn ( std::vector < double > & vec, int column )
00143 {
00144 int anynul;
00145 double nulval = 0;
00146 long nelements = vec.size();
00147 vector<double>::iterator it = vec.begin();
00148 double * ptr = &*it;
00149
00150 m_status = 0;
00151 fits_read_col_dbl( m_fptr, column+1, 1, 1, nelements, nulval,
00152 ptr, &anynul, &m_status);
00153
00154 return m_status;
00155 }
00156
00157 void
00158 FitsFile::
00159 fillShape ( std::vector < intptr_t > & shape, int column )
00160 {
00161 m_status = 0;
00162
00163 string tdimn ( "TDIM" );
00164 tdimn += String::convert ( column+1 );
00165 char value[80];
00166 fits_read_keyword ( m_fptr,
00167 const_cast < char * > ( tdimn.c_str() ),
00168 value,
00169 NULL, & m_status );
00170
00171 long rows = getNumberOfRows ();
00172 if ( value[0] != 0 ) {
00173 const string dims = value;
00174 unsigned int count = std::count ( dims.begin(), dims.end(), ',' );
00175 count += 1;
00176 int icount = 0;
00177 vector < long > naxes ( count );
00178 fits_read_tdim ( m_fptr, column+1, count,
00179 &icount, &naxes[0], & m_status );
00180 shape.resize ( count + 1 );
00181 shape [0] = rows;
00182 std::copy ( naxes.begin(), naxes.end(), shape.begin() + 1 );
00183 }
00184 else {
00185 shape.resize ( 1, rows );
00186 }
00187 }
00188
00189 int
00190 FitsFile::
00191 fillAxisSizes ( std::vector < long > & vec ) const
00192 {
00193 vec.clear ();
00194 int naxis = getImageDimensions ();
00195 vec.resize ( naxis );
00196
00197 vector < long > ::iterator first = vec.begin();
00198 long * ptr = & *first;
00199 m_status = 0;
00200 fits_get_img_size ( m_fptr, naxis, ptr, & m_status );
00201 assert ( m_status == 0 );
00202
00203 return m_status;
00204 }
00205
00206 void
00207 FitsFile::
00208 fillImageDeltas ( std::vector < double > & deltas ) const
00209 {
00210 deltas.clear ();
00211 int naxis = getImageDimensions ();
00212 deltas.reserve ( naxis );
00213
00214 char key [ FLEN_KEYWORD];
00215 char * keyroot = const_cast < char * > ("CDELT");
00216 for ( int i = 0; i < naxis; i++ ) {
00217 m_status = 0;
00218 fits_make_keyn ( keyroot, i+1, key, & m_status );
00219 bool yes = hasKey ( key );
00220
00221 if ( yes ) {
00222 double value = doubleValueForKey ( key );
00223 deltas.push_back ( value );
00224 }
00225 else {
00226 deltas.push_back ( 1.0 );
00227 }
00228 }
00229 }
00230
00231 void
00232 FitsFile::
00233 fillRefPixelIndices ( std::vector < int > & indices ) const
00234 {
00235 indices.clear ();
00236 int naxis = getImageDimensions ();
00237 indices.reserve ( naxis );
00238
00239 char key [ FLEN_KEYWORD];
00240 char * keyroot = const_cast < char * > ( "CRPIX" );
00241 for ( int i = 0; i < naxis; i++ ) {
00242 m_status = 0;
00243 fits_make_keyn ( keyroot, i+1, key, & m_status );
00244 bool yes = hasKey ( key );
00245
00246 if ( yes ) {
00247 int value = intValueForKey ( key );
00248 indices.push_back ( value );
00249 }
00250 else {
00251 indices.push_back ( 1 );
00252 }
00253 }
00254 }
00255
00256 void
00257 FitsFile::
00258 fillRefPixelValues ( std::vector < double > & values ) const
00259 {
00260 values.clear ();
00261 int naxis = getImageDimensions ();
00262 values.reserve ( naxis );
00263
00264 char key [ FLEN_KEYWORD];
00265 char * keyroot = const_cast < char * > ( "CRVAL" );
00266 for ( int i = 0; i < naxis; i++ ) {
00267 m_status = 0;
00268 fits_make_keyn ( keyroot, i+1, key, & m_status );
00269 bool yes = hasKey ( key );
00270
00271 if ( yes ) {
00272 double value = doubleValueForKey ( key );
00273 values.push_back ( value );
00274 }
00275 else {
00276 values.push_back ( 0. );
00277 }
00278 }
00279 }
00280
00281 bool
00282 FitsFile::
00283 isHammerAitoff () const
00284 {
00285 bool hammer = true;
00286
00287 char key [FLEN_KEYWORD];
00288 char * keyroot = const_cast < char * > ( "CTYPE" );
00289 for ( int i = 0; i < 2; i++ ) {
00290 fits_make_keyn ( keyroot, i+1, key, & m_status );
00291 bool yes = hasKey ( key );
00292
00293 if ( yes ) {
00294 string value = stringValueForKey ( key );
00295 if ( value.find ( "-AIT" ) == string::npos ) {
00296 hammer = false;
00297 }
00298 } else {
00299 hammer = false;
00300 }
00301 }
00302
00303 return hammer;
00304 }
00305
00306 int
00307 FitsFile::
00308 fillFromImage ( std::vector < double > & vec, unsigned int zplane )
00309 {
00310 vector < long > naxes;
00311 fillAxisSizes ( naxes );
00312
00313 int datatype = Double;
00314 long nelements = naxes[0] * naxes[1];
00315 double nulval = 0;
00316 vector < double >::iterator first = vec.begin();
00317 double * ptr = & *first;
00318 int anynul;
00319 m_status = 0;
00320 if ( naxes.size () == 2 ) {
00321 long fpixel[2] = { 1, 1 };
00322 fits_read_pix ( m_fptr, datatype, fpixel, nelements,
00323 & nulval, ptr, & anynul, & m_status );
00324 } else {
00325 long fpixel[] = { 1, 1, 1 };
00326 fpixel[2] = zplane + 1;
00327 fits_read_pix ( m_fptr, datatype, fpixel, nelements,
00328 & nulval, ptr, & anynul, & m_status );
00329 }
00330
00331 return m_status;
00332 }
00333
00334 int
00335 FitsFile::
00336 fillIntVectorFromColumn( std::vector < int > & vec, int column )
00337 {
00338 int anynul, status = 0;
00339 int nulval = 0;
00340 long nelements = vec.size();
00341 vector<int>::iterator it = vec.begin();
00342 int * ptr = new int [ nelements ];
00343
00344 m_status = 0;
00345 fits_read_col_int( m_fptr, column+1, 1, 1, nelements, nulval,
00346 ptr, &anynul, &status);
00347 copy ( ptr, ptr+nelements, it );
00348
00349 return status;
00350 }
00351
00352 void
00353 FitsFile::
00354 writeHDU ( long rows, int columns,
00355 const std::vector < std::string > & names,
00356 const std::vector < std::vector < int > > & shapes,
00357 const std::string & extname )
00358 {
00359 char ** types = new char * [ columns ];
00360 char ** tform = new char * [ columns ];
00361 char ** tunits = 0;
00362 char * m_extname;
00363
00364 vector < string > forms;
00365 for ( int i = 0; i < columns; i ++ ) {
00366 types[i]=const_cast<char*> (names[i].c_str());
00367 int size = 1;
00368 const vector < int > & shape = shapes [ i ];
00369 unsigned int rank = shape.size();
00370 if ( rank > 1 ) {
00371 for ( unsigned int j = 1; j < rank; j++ ) {
00372 int dim = shape [ j ];
00373 size *= dim;
00374 }
00375 }
00376 string form = String::convert ( size );
00377 form += "D";
00378 forms.push_back ( form );
00379
00380 tform[i] = strdup ( form.c_str() );
00381
00382
00383 }
00384
00385 m_extname = const_cast<char*> (extname.c_str());
00386
00387 fits_create_tbl( m_fptr, BINARY_TBL, 0, columns,
00388 types, tform, tunits, m_extname, &m_status);
00389
00390 for ( int i = 0; i < columns; i ++ ) {
00391 const vector < int > & shape = shapes [ i ];
00392 unsigned int rank = shape.size();
00393 if ( rank > 1 ) {
00394 vector < long > naxes;
00395 for ( unsigned int j = 1; j < rank; j++ ) {
00396 naxes.push_back ( shape [ j ] );
00397 }
00398 fits_write_tdim ( m_fptr, i+1,
00399 rank -1, &naxes [0], &m_status );
00400 }
00401 }
00402
00403 delete types;
00404 delete tform;
00405 }
00406
00407 void
00408 FitsFile::
00409 writeImageHDU ( long x, long y )
00410 {
00411 long naxes[2]={x,y};
00412 fits_create_img ( m_fptr, DOUBLE_IMG, 2, naxes, &m_status );
00413 }
00414
00415
00416 void
00417 FitsFile::
00418 writeColumn ( int col, const std::vector < double > & data )
00419 {
00420 LONGLONG nelements = data.size ();
00421 LONGLONG firstrow = 1;
00422 LONGLONG firstelem = 1;
00423 fits_write_col( m_fptr, TDOUBLE, col+1, firstrow, firstelem, nelements,
00424 const_cast < double * > ( &data[0] ), &m_status);
00425 }
00426
00427 void
00428 FitsFile::
00429 writePix ( long x, long y, const std::vector < double > & data )
00430 {
00431 long fpixel[2]={1,1};
00432 long nelements = x * y;
00433 fits_write_pix ( m_fptr, TDOUBLE, fpixel, nelements,
00434 const_cast < double * > ( &data[0] ), &m_status );
00435 }
00436
00437 void
00438 FitsFile::
00439 writeCloseFile ()
00440 {
00441 fits_close_file( m_fptr, &m_status );
00442 m_fptr = 0;
00443 }
00444
00445 bool
00446 FitsFile::
00447 pixCenter () const
00448 {
00449
00450 char * key = const_cast < char * > ( "PIXCENT" );
00451 bool yes = hasKey ( key );
00452 if ( yes )
00453 {
00454 if ( stringValueForKey ( key ) == "T" ) return true;
00455 }
00456
00457 return false;
00458
00459 }
00460
00461 void
00462 FitsFile::
00463 writeRefPixelValues ( double value1, double value2 )
00464 {
00465 char * key1 = const_cast < char * > ( "CRVAL1" );
00466 char * key2 = const_cast < char * > ( "CRVAL2" );
00467
00468 fits_write_key( m_fptr, TDOUBLE, key1, &value1, NULL, &m_status );
00469 fits_write_key( m_fptr, TDOUBLE, key2, &value2, NULL, &m_status );
00470 }