bes  Updated for version 3.20.6
fits_read_descriptors.cc
1 // fits_read_descriptors.cc
2 
3 // This file is part of fits_handler, a data handler for the OPeNDAP data
4 // server.
5 
6 // Copyright (c) 2004,2005 University Corporation for Atmospheric Research
7 // Author: Patrick West <pwest@ucar.edu> and Jose Garcia <jgarcia@ucar.edu>
8 //
9 // This library is free software; you can redistribute it and/or
10 // modify it under the terms of the GNU Lesser General Public
11 // License as published by the Free Software Foundation; either
12 // version 2.1 of the License, or (at your option) any later version.
13 //
14 // This library is distributed in the hope that it will be useful,
15 // but WITHOUT ANY WARRANTY; without even the implied warranty of
16 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 // Lesser General Public License for more details.
18 //
19 // You should have received a copy of the GNU Lesser General Public
20 // License along with this library; if not, write to the Free Software
21 // Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
22 //
23 // You can contact University Corporation for Atmospheric Research at
24 // 3080 Center Green Drive, Boulder, CO 80301
25 
26 // (c) COPYRIGHT University Corporation for Atmospheric Research 2004-2005
27 // Please read the full copyright statement in the file COPYRIGHT_UCAR.
28 //
29 // Authors:
30 // pwest Patrick West <pwest@ucar.edu>
31 // jgarcia Jose Garcia <jgarcia@ucar.edu>
32 
33 #include "config.h"
34 
35 #include <sstream>
36 #include <memory>
37 
38 #include <DDS.h>
39 #include <Structure.h>
40 #include <Str.h>
41 #include <Array.h>
42 #include <Byte.h>
43 #include <Int16.h>
44 #include <Int32.h>
45 #include <Float32.h>
46 #include <Float64.h>
47 #include <Error.h>
48 #include <mime_util.h>
49 
50 #include <BESDebug.h>
51 
52 #include <fitsio.h>
53 
54 #include "fits_read_descriptors.h"
55 #include "BESAutoPtr.h"
56 
57 using namespace std;
58 using namespace libdap;
59 
60 void fits_handler::process_status(int status, string &error)
61 {
62  if (status) {
63  fits_report_error(stderr, status);
64  }
65  char error_description[30] = "";
66  fits_get_errstatus(status, error_description);
67  error = error_description;
68  return;
69 }
70 
71 #if 0
72 char *
73 fits_handler::ltoa(long val, /* value to be converted */
74 char *buf, /* output string */
75 int base) /* conversion base */
76 {
77  ldiv_t r; /* result of val / base */
78 
79  if (base > 36 || base < 2) /* no conversion if wrong base */
80  {
81  *buf = '\0';
82  return buf;
83  }
84  if (val < 0)
85  *buf++ = '-';
86  r = ldiv(labs(val), base);
87 
88  /* output digits of val/base first */
89 
90  if (r.quot > 0)
91  buf = ltoa(r.quot, buf, base);
92  /* output last digit */
93 
94  *buf++ = "0123456789abcdefghijklmnopqrstuvwxyz"[(int) r.rem];
95  *buf = '\0';
96  return buf;
97 }
98 #endif
99 
100 bool fits_handler::fits_read_descriptors(DDS &dds, const string &filename, string &error)
101 {
102  fitsfile *fptr;
103  int status = 0;
104  if (fits_open_file(&fptr, filename.c_str(), READONLY, &status)) {
105  error = "Can not open fits file ";
106  error += filename;
107  return false;
108  }
109 
110  dds.set_dataset_name(name_path(filename));
111 
112  int hdutype;
113  for (int ii = 1; !(fits_movabs_hdu(fptr, ii, &hdutype, &status)); ii++) {
114  ostringstream hdu;
115  hdu << "HDU_" << ii;
116  switch (hdutype) {
117  case IMAGE_HDU:
118  // 'hdu' holds the name of the Structure; hdu + _IMAGE is the name of
119  // the variable within that structure that holds the image data.
120  status = process_hdu_image(fptr, dds, hdu.str(), hdu.str() + "_IMAGE");
121  process_status(status, error);
122  break;
123  case ASCII_TBL:
124  status = process_hdu_ascii_table(fptr, dds, hdu.str(), hdu.str() + "_ASCII");
125  process_status(status, error);
126  break;
127  case BINARY_TBL:
128  // FIXME This should return an error!
129  // process_hdu_binary_table does/did nothing but return a non-error status code.
130  // NB: I rewrote the call, commented out, using the two new args that eliminate
131  // the global string variables. jhrg 6/14/13
132  status = 0; // process_hdu_binary_table(fptr, dds, hdu.str(), hdu.str() + "_BINARY");
133  process_status(status, error);
134  break;
135  default:
136  // cerr << "Invalid HDU type in file " << filename << endl;
137  process_status(1, error);
138  break;
139  }
140  }
141 
142  if (status == END_OF_FILE) // status values are defined in fitsioc.h
143  status = 0; // got the expected EOF error; reset = 0
144  else {
145  process_status(status, error);
146  fits_close_file(fptr, &status);
147  return false;
148  }
149  if (fits_close_file(fptr, &status)) {
150  process_status(status, error);
151  return false;
152  }
153  return true;
154 }
155 
156 // T = Byte --> U = dods_byte, fits_types = TBYTE
157 template<class T, class U>
158 static int read_hdu_image_data(int fits_type, fitsfile *fptr, const string &varname, const string &datasetname, int number_axes, const vector<long> &naxes, Structure *container)
159 {
160  unique_ptr<T> in(new T(varname, datasetname));
161  unique_ptr<Array> arr(new Array(varname, datasetname, in.get()));
162  long npixels = 1;
163  for (register int w = 0; w < number_axes; w++) {
164  ostringstream oss;
165  oss << "NAXIS" << w;
166  arr->append_dim(naxes[w], oss.str());
167 
168  npixels = npixels * naxes[w];
169  }
170  dods_byte nullval = 0;
171  vector<U> buffer(npixels);
172  int anynull, status = 0;
173  // the original code always read the whole array; I dropped passing in fpixel and replaced it with a '1'.
174  // jhrg 6/14/13
175  fits_read_img(fptr, fits_type, 1 /*first pixel*/, npixels, &nullval, &buffer[0], &anynull, &status);
176  arr->set_value(&buffer[0], npixels);
177  container->add_var(arr.get());
178 
179  return status;
180 }
181 
182 int fits_handler::process_hdu_image(fitsfile *fptr, DDS &dds, const string &hdu, const string &str)
183 {
184  string datasetname = dds.get_dataset_name();
185 
186  unique_ptr<Structure> container(new Structure(hdu, datasetname));
187 
188  int status = 0;
189  int nkeys, keypos;
190  if (fits_get_hdrpos(fptr, &nkeys, &keypos, &status))
191  return status;
192 
193  for (int jj = 1; jj <= nkeys; jj++) {
194  char name[FLEN_KEYWORD];
195  char value[FLEN_VALUE];
196  char comment[FLEN_COMMENT];
197  if (fits_read_keyn(fptr, jj, name, value, comment, &status)) return status;
198 
199  ostringstream keya;
200  keya << "key_" << jj;
201  unique_ptr<Structure> st(new Structure(keya.str(), datasetname));
202 
203  unique_ptr<Str> s1(new Str("name", datasetname));
204  string ppp = name;
205  s1->set_value(ppp);
206 
207  unique_ptr<Str> s2(new Str("value", datasetname));
208  ppp = value;
209  s2->set_value(ppp);
210 
211  unique_ptr<Str> s3(new Str("comment", datasetname));
212  ppp = comment;
213  s3->set_value(ppp);
214 
215  st->add_var(s1.get());
216  st->add_var(s2.get());
217  st->add_var(s3.get());
218  container->add_var(st.get());
219  }
220 
221  char value[FLEN_VALUE];
222  if (fits_read_keyword(fptr, "BITPIX", value, NULL, &status))
223  return status; // status is not 0, so something is wrong and we get until here...
224  int bitpix = atoi(value);
225 
226  if (fits_read_keyword(fptr, "NAXIS", value, NULL, &status))
227  return status;
228  int number_axes = atoi(value);
229 
230  vector<long> naxes(number_axes);
231  int nfound;
232  if (fits_read_keys_lng(fptr, "NAXIS", 1, number_axes, &naxes[0], &nfound, &status))
233  return status; // status is not 0, so something is wrong
234 
235  switch (bitpix) {
236  case BYTE_IMG:
237  status = read_hdu_image_data<Byte,dods_byte>(TBYTE, fptr, str, datasetname, number_axes, naxes, container.get());
238  break;
239 
240  case SHORT_IMG:
241  status = read_hdu_image_data<Int16,dods_int16>(TSHORT, fptr, str, datasetname, number_axes, naxes, container.get());
242  break;
243 
244  case LONG_IMG:
245  // Here is/was the bug that started me down the whole rabbit hole of reworking this code. Passing TLONG
246  // instead of TINT seems to make cfitsio 3270 to 3340 (and others?) treat the data as if it is an
247  // array of 64-bit integers. This is odd since TLONG is a 32-bit int type #define'd in cfitsio.h.
248  // In valgrind, using TLONG returns the following stack trace from valgrind:
249  //
250  // ==22851== Invalid write of size 8
251  // ==22851== at 0x3186835: fffi4i4 (getcolj.c:1158)
252  // ==22851== by 0x31899A2: ffgclj (getcolj.c:779)
253  // ==22851== by 0x3186372: ffgpvj (getcolj.c:56)
254  // ==22851== by 0x3178719: ffgpv (getcol.c:637)
255  // ==22851== by 0x3122127: fits_handler::process_hdu_image(fitsfile*, libdap::DDS&, std::string const&, std::string const&) (fits_read_descriptors.cc:169)
256  //
257  // valgrind --dsymutil=yes besstandalone -c bes.conf -i fits/dpm.dds.bescmd
258  //
259  // If TLONG is used and the array of dods_int32 is allocated as 2 * npixels, the code does not trigger
260  // the invalid write message.
261  //
262  // jhrg 6/14/13
263  status = read_hdu_image_data<Int32,dods_int32>(TINT/*TLONG*/, fptr, str, datasetname, number_axes, naxes, container.get());
264  break;
265 
266  case FLOAT_IMG:
267  status = read_hdu_image_data<Float32,dods_float32>(TFLOAT, fptr, str, datasetname, number_axes, naxes, container.get());
268  break;
269 
270  case DOUBLE_IMG:
271  status = read_hdu_image_data<Float64,dods_float64>(TDOUBLE, fptr, str, datasetname, number_axes, naxes, container.get());
272  break;
273 
274  default:
275  status = 1;
276  break;
277  }
278 
279  dds.add_var(container.get());
280  return status;
281 }
282 
283 // I made minimal changes to this code below; mostly using auto_ptr where applicable and fixing
284 // the bug with cfitsio where TLONG seems to be treated as a 64-but integer.
285 
286 int fits_handler::process_hdu_ascii_table(fitsfile *fptr, DDS &dds, const string &hdu, const string &str)
287 {
288  string datasetname = dds.get_dataset_name();
289  unique_ptr<Structure> container(new Structure(hdu, datasetname));
290  int status = 0;
291  int nfound, anynull;
292  int ncols;
293  long nrows;
294  int nkeys, keypos;
295  char name[FLEN_KEYWORD];
296  char value[FLEN_VALUE];
297  char comment[FLEN_COMMENT];
298  anynull = 0;
299 
300  if (fits_get_hdrpos(fptr, &nkeys, &keypos, &status))
301  return status;
302 
303  for (int jj = 1; jj <= nkeys; jj++) {
304  if (fits_read_keyn(fptr, jj, name, value, comment, &status))
305  return status;
306  ostringstream oss;
307  oss << "key_" << jj;
308  unique_ptr<Structure> st(new Structure(oss.str(), datasetname));
309  unique_ptr<Str> s1(new Str("name", datasetname));
310  string ppp = name;
311  s1->set_value(ppp);
312  unique_ptr<Str> s2(new Str("value", datasetname));
313  ppp = value;
314  s2->set_value(ppp);
315  unique_ptr<Str> s3(new Str("comment", datasetname));
316  ppp = comment;
317  s3->set_value(ppp);
318  st->add_var(s1.get());
319  st->add_var(s2.get());
320  st->add_var(s3.get());
321  container->add_var(st.get());
322  }
323 
324  fits_get_num_rows(fptr, &nrows, &status);
325  fits_get_num_cols(fptr, &ncols, &status);
326 
327  // I am sure this is one of the most obscure piece of code you have ever seen
328  // First I get an auto pointer object to hold securely an array of auto pointer
329  // objects holding securely pointers to char...
330  BESAutoPtr<BESAutoPtr<char> > ttype_autoptr(new BESAutoPtr<char> [ncols], true);
331 
332  // Now I set each one of my auto pointer objects holding pointers to char
333  // to hold the address of a dynamically allocated piece of memory (array of chars)...
334  for (int j = 0; j < ncols; j++) {
335  (ttype_autoptr.get())[j].reset();
336  (ttype_autoptr.get())[j].set(new char[FLEN_VALUE], true);
337  }
338 
339  // Now I align my pointers to char inside each BESAutoPtr <char> object
340  // inside my BESAutoPtr< BESAutoPtr <char> > object to a char** pointer.
341 
342  // Step 1:
343  // I get a insecure pointer an get some memory on it...
344  char **ttype = new char*[ncols];
345 
346  // Step 2;
347  // Lets secure the memory area pointed by ttype inside
348  // a BESAutoPtr<char*> object, Lets not forget ttype is a vector
349  BESAutoPtr<char*> secure_ttype(ttype, true);
350 
351  // Step 3:
352  // OK here we are, now we have ncols beautifully aligned pointers to char
353  // let the pointer inside ttype point to the same address as the secure ones...
354  for (int j = 0; j < ncols; j++)
355  ttype[j] = (ttype_autoptr.get())[j].get();
356 
357  // Step 4:
358  // Now we read the data!
359  if (fits_read_keys_str(fptr, "TTYPE", 1, ncols, ttype, &nfound, &status))
360  return status;
361 
362  // wasn't that fun ? :)
363 
364 
365  unique_ptr<Structure> table(new Structure(str, datasetname));
366 
367  for (int h = 0; h < ncols; h++) {
368  int typecode;
369  long repeat, width;
370  fits_get_coltype(fptr, h + 1, &typecode, &repeat, &width, &status);
371 
372  switch (typecode) {
373  case TSTRING: {
374  int p;
375  unique_ptr<Str> in(new Str(ttype[h], datasetname));
376  unique_ptr<Array> arr(new Array(ttype[h], datasetname, in.get()));
377  arr->append_dim(nrows);
378  char strnull[10] = "";
379  char **name = new char*[nrows];
380  // secure the pointer for exceptions, remember is an array so second parameter is true
381  BESAutoPtr<char *> sa1(name, true);
382  // get an auto pointer (sa3) object to hold securely an array of auto pointers to char
383  BESAutoPtr<BESAutoPtr<char> > sa3(new BESAutoPtr<char> [nrows], true);
384  for (p = 0; p < nrows; p++) {
385  // get memory
386  name[p] = new char[width + 1];
387  // reset auto pointer
388  (sa3.get())[p].reset();
389  // storage safely the area in the heap pointed by name[p] in an auto pointer
390  (sa3.get())[p].set(name[p], true);
391  }
392  fits_read_col(fptr, typecode, h + 1, 1, 1, nrows, strnull, name, &anynull, &status);
393 
394  string *strings = new string[nrows];
395  // secure the pointer for exceptions, remember is an array
396  BESAutoPtr<string> sa2(strings, true);
397 
398  for (int p = 0; p < nrows; p++)
399  strings[p] = name[p];
400  arr->set_value(strings, nrows);
401 
402  table->add_var(arr.get());
403  }
404  break;
405  case TSHORT: {
406  BESAutoPtr<Int16> in(new Int16(ttype[h], datasetname));
407  BESAutoPtr<Array> arr(new Array(ttype[h], datasetname, in.get()));
408  arr->append_dim(nrows);
409  dods_int16 nullval = 0;
410  BESAutoPtr<dods_int16> buffer(new dods_int16[nrows], true);
411  fits_read_col(fptr, typecode, h + 1, 1, 1, nrows, &nullval, buffer.get(), &anynull, &status);
412  arr->set_value(buffer.get(), nrows);
413  table->add_var(arr.get());
414  }
415  break;
416  case TLONG: {
417  BESAutoPtr<Int32> in(new Int32(ttype[h], datasetname));
418  BESAutoPtr<Array> arr(new Array(ttype[h], datasetname, in.get()));
419  arr->append_dim(nrows);
420  dods_int32 nullval = 0;
421  BESAutoPtr<dods_int32> buffer(new dods_int32[nrows], true);
422  fits_read_col(fptr, TINT/*typecode*/, h + 1, 1, 1, nrows, &nullval, buffer.get(), &anynull, &status);
423  arr->set_value(buffer.get(), nrows);
424  table->add_var(arr.get());
425  }
426  break;
427  case TFLOAT: {
428  BESAutoPtr<Float32> in(new Float32(ttype[h], datasetname));
429  BESAutoPtr<Array> arr(new Array(ttype[h], datasetname, in.get()));
430  arr->append_dim(nrows);
431  dods_float32 nullval = 0;
432  BESAutoPtr<dods_float32> buffer(new dods_float32[nrows], true);
433  fits_read_col(fptr, typecode, h + 1, 1, 1, nrows, &nullval, buffer.get(), &anynull, &status);
434  arr->set_value(buffer.get(), nrows);
435  table->add_var(arr.get());
436  }
437  break;
438  case TDOUBLE: {
439  BESAutoPtr<Float64> in(new Float64(ttype[h], datasetname));
440  BESAutoPtr<Array> arr(new Array(ttype[h], datasetname, in.get()));
441  arr->append_dim(nrows);
442  dods_float64 nullval = 0;
443  BESAutoPtr<dods_float64> buffer(new dods_float64[nrows], true);
444  fits_read_col(fptr, typecode, h + 1, 1, 1, nrows, &nullval, buffer.get(), &anynull, &status);
445  arr->set_value(buffer.get(), nrows);
446  table->add_var(arr.get());
447  }
448  break;
449  }
450  }
451  container->add_var(table.get());
452  dds.add_var(container.get());
453  return status;
454 }
455 
456 #if 0
457 int fits_handler::process_hdu_binary_table(fitsfile *, DDS &)
458 {
459  return 0;
460 }
461 #endif
BESAutoPtr
Definition: BESAutoPtr.h:37
libdap
Definition: BESDapFunctionResponseCache.h:35