mfitsio.c Source File

00001 
00030 #include <string.h>
00031 #include <stdio.h>
00032 #include <fitsio.h>
00033 #include <malloc.h>
00034 #include <unistd.h>
00035 
00037 #include "mex.h"
00038 #include "matrix.h"
00039 
00041 #include "mfitsio.h"
00042 
00043 /*************************************************************************
00044  *  Allocation and deallocation functions for mfitsio data structures. *
00045  *************************************************************************/
00046 
00053 void
00054 mfitsio_free_record (mfitsio_record * record)
00055 {
00056   MFITSIO_FREE (record->key);
00057   MFITSIO_FREE (record->value);
00058 }
00059 
00066 void
00067 mfitsio_free_header (mfitsio_header * header)
00068 {
00069   int i;
00070   for (i = 0; i < header->length; i++)
00071   {
00072     mfitsio_free_record (header->records[i]);
00073   }
00074   MFITSIO_FREE (header->records);
00075 }
00076 
00083 void
00084 mfitsio_free_info (mfitsio_info * info)
00085 {
00086   MFITSIO_FREE (info->naxes);
00087   MFITSIO_FREE (info);
00088 }
00089 
00103 mfitsio_record *
00104 mfitsio_parse_card (const char *cardtext)
00105 {
00106   int i = 0, j = 0, k = 0, len_card = 0, len_key = 0, len_value = 0;
00107   char value = 0;
00108   mfitsio_record *record = 0;
00109 
00111   if (!mfitsio_ignore_card (cardtext))
00112   {
00114     record = (mfitsio_record *) MFITSIO_MALLOC (sizeof (mfitsio_record));
00115     record->key = 0;
00116     record->value = 0;
00117 
00119     len_card = strlen (cardtext);
00120 
00123     for (i = 0; i < len_card; i++)
00124     {
00126       if (cardtext[i] == ' ')
00127       {
00128         break;
00129       }
00133       else if (cardtext[i] == '=')
00134       {
00135         k = 1;
00136         break;
00137       }
00138     }
00139 
00141     len_key = i + 1;
00142     record->key = (char *) MFITSIO_MALLOC (len_key * sizeof (char));
00143 
00145     memcpy (record->key, cardtext, i);
00146 
00148     record->key[i] = 0;
00149 
00151     j = i + 1;
00152 
00155     if (k == 0)
00156     {
00157       for (; j < len_card; j++)
00158       {
00159         if (cardtext[j] == '=')
00160         {
00161           j++;
00162           break;
00163         }
00164       }
00165     }
00166 
00168     for (; j < len_card; j++)
00169     {
00170       if (cardtext[j] != ' ')
00171       {
00172         break;
00173       }
00174     }
00175 
00177     len_value = (len_card - j + 1);
00178 
00180     record->value = (char *) MFITSIO_MALLOC (len_value * sizeof (char));
00181     memcpy (record->value, cardtext + j, len_value - 1);
00182 
00184     record->value[len_value - 1] = 0;
00185 
00191     for (j = len_value - 1; j >= 0; j--)
00192     {
00193       value = record->value[j];
00194       if (value == '\'')
00195       {
00196         break;
00197       }
00198       else if (value == '/' || value == ' ')
00199       {
00200         record->value[j] = '\0';
00201       };
00202     }
00203   }
00204 
00208   return record;
00209 }
00210 
00222 mfitsio_header *
00223 mfitsio_read_header (const char *filename)
00224 {
00225   fitsfile *fptr = 0;
00226   mfitsio_header *header = 0, *new_header = 0;
00227   char card[FLEN_CARD];
00228   int status = 0, nkeys = 0, jkeys = 0, i = 0, j = 0;
00229   mfitsio_record *record = 0;
00230 
00232   fits_open_file (&fptr, filename, MFITSIO_READONLY, &status);
00233   if (status)
00234   {
00235     MFITSIO_PRINTF ("1  An error occurred: %d\n", status);
00236     MFITSIO_ERR ("Unable to open FITS file.");
00237   }
00238 
00240   fits_get_hdrspace (fptr, &nkeys, NULL, &status);
00241   if (status)
00242   {
00243     MFITSIO_PRINTF ("2  An error occurred: %d\n", status);
00244     MFITSIO_ERR ("Could not access header.");
00245   }
00246 
00248   header = (mfitsio_header *) MFITSIO_MALLOC (sizeof (mfitsio_header));
00249 
00251   header->length = nkeys;
00252 
00254   header->records =
00255     (mfitsio_record **) MFITSIO_MALLOC (nkeys * sizeof (mfitsio_record *));
00256 
00259   for (i = 1; i <= nkeys; i++)
00260   {
00261     fits_read_record (fptr, i, card, &status);
00262     if (status)
00263     {
00264       MFITSIO_PRINTF ("3  An error occurred: %d\n", status);
00265       MFITSIO_ERR ("Could not access header.");
00266     }
00267     record = mfitsio_parse_card (card);
00268     header->records[i - 1] = record;
00269 
00271     if (record != 0)
00272     {
00273       jkeys++;
00274     }
00275   }
00276 
00280   if (jkeys != nkeys)
00281   {
00283     new_header =
00284       (mfitsio_header *) MFITSIO_MALLOC (sizeof (mfitsio_header));
00285     new_header->records =
00286       (mfitsio_record **) MFITSIO_MALLOC (jkeys *
00287                                             sizeof (mfitsio_record *));
00288 
00290     new_header->length = jkeys;
00291     j = 0;
00292     for (i = 0; i < nkeys; i++)
00293     {
00294       if (header->records[i] != 0)
00295       {
00296         new_header->records[j++] = header->records[i];
00297       }
00298     }
00300     MFITSIO_FREE (header->records);
00301     MFITSIO_FREE (header);
00302 
00304     header = new_header;
00305   }
00306 
00308   fits_close_file (fptr, &status);
00309   if (status)
00310   {
00311     MFITSIO_PRINTF ("4  An error occurred: %d\n", status);
00312     MFITSIO_ERR ("Could not close file.");
00313   }
00314 
00315   return (status != 0) ? 0 : header;
00316 }
00317 
00327 mxArray *
00328 mfitsio_adapt_frecord (const mfitsio_record * record)
00329 {
00330   int len = 0, i, j;
00331   char *temp = 0;
00332   char f = 0;
00333 
00335   mxArray *retval = 0;
00336   mxArray *iargs[] = { 0 };
00337   mxArray *oargs[] = { 0 };
00338 
00340   len = strlen (record->value);
00341 
00343   if (len)
00344   {
00345 
00347     f = record->value[0];
00348 
00351     if (f == '\'')
00352     {
00354       temp = (char *) MFITSIO_MALLOC ((len - 1) * sizeof (char));
00355 
00358       memcpy (temp, record->value + 1, len - 2);
00359 
00361       temp[len - 2] = 0;
00362 
00365       i = -1;
00366       for (j = strlen (temp) - 1; j >= 0; j--)
00367       {
00368         if (temp[j] != ' ')
00369         {
00370           i = j + 1;
00371           break;
00372         }
00373       }
00374       if (i != -1)
00375       {
00376         temp[i] = '\0';
00377       }
00378 
00380       retval = mxCreateString (temp);
00381     }
00383     else if (f == 'T' || f == 'F' || f == 't' || f == 'f')
00384     {
00385       retval = mxCreateLogicalScalar (f == 'T' || f == 't');
00386     }
00388     else if ((f >= '0' && f <= '9') || f == '-' || f == '.')
00389     {
00392       if (strchr (record->value, '.'))
00393       {
00394         retval = mxCreateDoubleScalar (atof (record->value));
00395       }
00397       else
00398       {
00400         retval = mxCreateDoubleScalar (atof (record->value));
00401         iargs[0] = retval;
00402 
00404         mexCallMATLAB (1, oargs, 1, iargs, "int32");
00405         retval = oargs[0];
00406       }
00407     }
00409     else
00410     {
00411       retval = mxCreateString (strdup (record->value));
00412       i = -1;
00413       for (j = strlen (record->value) - 1; j >= 0; j--)
00414       {
00415         if (record->value[j] != ' ')
00416         {
00417           i = j + 1;
00418           break;
00419         }
00420       }
00421       if (i != -1)
00422       {
00423         record->value[i] = '\0';
00424       }
00425     }
00426   }
00428   else
00429   {
00430     retval = mxCreateDoubleScalar (0);
00431   }
00432   return retval;
00433 }
00434 
00444 mxArray *
00445 mfitsio_adapt_fheader (const mfitsio_header * header)
00446 {
00447   const char **fnames = 0;
00448   int buflen, status = 0, i;
00449   char *temp = 0;
00450   mxArray *retval = 0;
00451 
00452   if (header == 0)
00453   {
00454     MFITSIO_ERR ("33 Header unreadable.");
00455   }
00457   fnames = MFITSIO_CALLOC (header->length, sizeof (*fnames));
00458 
00461   for (i = 0; i < header->length; i++)
00462   {
00464     buflen = strlen (header->records[i]->key);
00465 
00467     temp = (char *) MFITSIO_MALLOC (sizeof (char) * buflen + 1);
00468 
00470     memcpy (temp, header->records[i]->key, buflen);
00471 
00473     temp[buflen] = 0;
00474 
00476     fnames[i] = temp;
00477   }
00478 
00480   retval = mxCreateStructMatrix (1, 1, header->length, fnames);
00481 
00486   for (i = 0; i < header->length; i++)
00487   {
00488     mxSetFieldByNumber (retval, 0, i,
00489                         mfitsio_adapt_frecord (header->records[i]));
00490   }
00491   return retval;
00492 }
00493 
00502 mfitsio_info *
00503 mfitsio_read_info (const char * filename)
00504 {
00505   mfitsio_info *retval = 0;
00506   fitsfile *fptr = 0;
00507   long *naxes;
00508   int status = 0, i;
00509 
00511   fits_open_file (&fptr, filename, MFITSIO_READONLY, &status);
00512   if (status)
00513   {
00514     MFITSIO_PRINTF ("5  An error occurred: %d\n", status);
00515     MFITSIO_ERR ("Could not open file for reading.");
00516   }
00517 
00519   retval = (mfitsio_info *) MFITSIO_MALLOC (sizeof (mfitsio_info));
00520 
00522   fits_get_img_type(fptr, &(retval->bitpix), &status);
00523   if (status)
00524   {
00525     MFITSIO_PRINTF ("6  An error occurred: %d\n", status);
00526     MFITSIO_ERR ("Could not retrieve BITPIX keyword.");
00527   }
00528 
00530   fits_get_img_dim(fptr, &(retval->naxis), &status);
00531   if (status)
00532   {
00533     MFITSIO_PRINTF ("7  An error occurred: %d\n", status);
00534     MFITSIO_ERR ("Could not retrieve NAXIS, total number of dimensions.");
00535   }
00536 
00538   naxes = MFITSIO_MALLOC (sizeof (long) * retval->naxis);
00539 
00541   fits_get_img_size(fptr, retval->naxis, naxes, &status);
00542   if (status)
00543   {
00544     MFITSIO_PRINTF ("8  An error occurred: %d\n", status);
00545     MFITSIO_ERR ("Could not open file for reading.");
00546   }
00547   retval->naxes = (int *) MFITSIO_MALLOC (retval->naxis * sizeof (int));
00548 
00550   for (i = 1; i <= retval->naxis; i++)
00551   {
00552     retval->naxes[i - 1] = (int) naxes[i - 1];
00553   }
00554 
00556   MFITSIO_FREE ( naxes );
00557 
00559   return retval;
00560 }
00561 
00571 mxArray *
00572 mfitsio_read_image (const char *filename, const mfitsio_info * info) {
00573   mxArray *image = 0;
00574   long *fpixels = 0, *lpixels = 0, *lnaxes = 0, *inc = 0;
00575   int i = 0;
00576   fpixels = (long *) MFITSIO_MALLOC (sizeof (long) * info->naxis);
00577   lpixels = (long *) MFITSIO_MALLOC (sizeof (long) * info->naxis);
00578   lnaxes = (long *) MFITSIO_MALLOC (sizeof (long) * info->naxis);
00579   inc = (long *) MFITSIO_MALLOC (sizeof (long) * info->naxis);
00581   for (i = 0; i < info->naxis; i++)
00582     {
00583       lpixels[i] = info->naxes[i];
00584       fpixels[i] = 1;
00585       inc[i] = 1;
00586       lnaxes[i] = info->naxes[i];
00587     }
00588   image = mfitsio_read_image_impl(filename, info, fpixels,
00589                                   lpixels, lnaxes, inc);
00590   MFITSIO_FREE (fpixels);
00591   MFITSIO_FREE (lpixels);
00592   MFITSIO_FREE (inc);
00593   MFITSIO_FREE (lnaxes);
00594   return image;
00595 }
00596 
00606 mxArray *
00607 mfitsio_read_image_impl (const char *filename, const mfitsio_info * info,
00608                          long *fpixels, long *lpixels,
00609                          long *lnaxes, long *inc)
00610 {
00611   fitsfile *fptr = 0;
00612   mxArray *img = 0;
00613   double *pr = 0;
00614   int status = 0, i = 0;
00615   int *diff = 0;
00616 
00617   diff = mfitsio_get_region_size(fpixels, lpixels, info->naxis);
00621   fits_open_file (&fptr, filename, MFITSIO_READONLY, &status);
00622   if (status)
00623   {
00624     MFITSIO_PRINTF ("9  An error occurred: %d\n", status);
00625     MFITSIO_ERR ("Could not open file for reading.");
00626   }
00627 
00629   img = mfitsio_read_image_execute(info, diff, fptr, 0, info->naxis, lnaxes,
00630                                    fpixels, lpixels, inc, 0, &status);
00631 
00632   MFITSIO_FREE (diff);
00633 
00635   fits_close_file (fptr, &status);
00636   if (status)
00637   {
00638     MFITSIO_PRINTF ("12 An error occurred: %d\n", status);
00639     MFITSIO_ERR ("Could not close file.");
00640   }
00641 
00643   return img;
00644 }
00645 
00646 mxArray *mfitsio_read_image_execute(const mfitsio_info *info, int *diff,
00647                                     fitsfile *fptr,
00648                                     int group, int naxis, long *lnaxes,
00649                                     long *fpixel, long *lpixel,
00650                                     long *inc, int *anynul, int *status) {
00651   mxArray *img;
00652   void *data;
00653   switch (info->bitpix)
00654     {
00655     case BYTE_IMG:
00656       img = mxCreateNumericArray (info->naxis, diff, mxUINT8_CLASS, mxREAL);
00657       data = mxGetData(img);
00658       fits_read_subset_byt (fptr, group, info->naxis, lnaxes, fpixel, lpixel,
00659                             inc, 0, (unsigned char*)data, anynul, status);
00660       break;
00661     case SHORT_IMG:
00662       img = mxCreateNumericArray (info->naxis, diff, mxINT16_CLASS, mxREAL);
00663       data = mxGetData(img);
00664       fits_read_subset_sht (fptr, 0, info->naxis, lnaxes, fpixel, lpixel,
00665                             inc, 0, (short*)data, 0, status);
00666       break;
00667     case LONG_IMG:
00668       img = mxCreateNumericArray (info->naxis, diff, mxINT32_CLASS, mxREAL);
00669       data = mxGetData(img);
00670       fits_read_subset_int (fptr, 0, info->naxis, lnaxes, fpixel, lpixel,
00671                             inc, 0, (int*)data, 0, status);
00672       break;
00673     case FLOAT_IMG:
00674       img = mxCreateNumericArray (info->naxis, diff, mxSINGLE_CLASS, mxREAL);
00675       data = mxGetData(img);
00676       fits_read_subset_flt (fptr, 0, info->naxis, lnaxes, fpixel, lpixel,
00677                             inc, 0, (float*)data, 0, status);
00678       break;
00679     case DOUBLE_IMG:
00680       img = mxCreateNumericArray (info->naxis, diff, mxDOUBLE_CLASS, mxREAL);
00681       data = mxGetData(img);
00682       fits_read_subset_dbl (fptr, 0, info->naxis, lnaxes, fpixel, lpixel,
00683                             inc, 0, (double*)data, 0, status);
00684       break;
00685     default:
00686       MFITSIO_PRINTF ("52 Unsupported input BITPIX: %d\n", info->bitpix);
00687       MFITSIO_ERR ("Cannot continue.");
00688     }  
00689   return img;
00690 }
00691 
00692 
00693 
00703 void
00704 mfitsio_write_image (const char *filename, const mxArray * header,
00705                       mxArray * img)
00706 {
00707 
00708   mfitsio_info *info = 0;
00709   long *fpixels = 0, *lpixels = 0, *lnaxes = 0, *inc = 0;
00710   int i;
00712   info = mfitsio_calc_info (img);
00713 
00715   fpixels = (long *) MFITSIO_MALLOC (sizeof (long) * info->naxis);
00716   lpixels = (long *) MFITSIO_MALLOC (sizeof (long) * info->naxis);
00717   lnaxes = (long *) MFITSIO_MALLOC (sizeof (long) * info->naxis);
00718   inc = (long *) MFITSIO_MALLOC (sizeof (long) * info->naxis);
00719   
00721   for (i = 0; i < info->naxis; i++)
00722     {
00723       lpixels[i] = info->naxes[i];
00724       fpixels[i] = 1;
00725       inc[i] = 1;
00726       lnaxes[i] = info->naxes[i];
00727     }
00728 
00729   mfitsio_write_image_impl(filename, header, img, fpixels,
00730                            lpixels, lnaxes, inc, info);
00731   MFITSIO_FREE (fpixels);
00732   MFITSIO_FREE (lpixels);
00733   MFITSIO_FREE (inc);
00734   MFITSIO_FREE (lnaxes);
00735   mfitsio_free_info(info);
00736 }
00737 
00750 void
00751 mfitsio_write_image_impl (const char *filename, const mxArray * header,
00752                           mxArray * img, long *fpixels, long *lpixels,
00753                           long *lnaxes, long *inc, mfitsio_info *info)
00754 {
00755 
00756   fitsfile *fptr = 0;
00757   mxArray *newImg = 0;
00758   int status = 0, nkeys, i, j;
00759   double *pr = 0;
00760 
00762   mfitsio_write_header (filename, header);
00763 
00765   fits_open_file (&fptr, filename, MFITSIO_READWRITE, &status);
00766   if (status)
00767   {
00768     MFITSIO_PRINTF ("13 An error occurred: %d\n", status);
00769     MFITSIO_ERR ("Could not open file for writing.");
00770   }
00771 
00773   newImg = img;
00791   if (info->check == 0) {
00792     fits_resize_img (fptr, info->bitpix, info->naxis, lnaxes, &status);
00793   }
00794   if (status)
00795   {
00796     MFITSIO_PRINTF ("14 An error occurred: %d\n", status);
00797     MFITSIO_ERR ("Could not resize FITS file image.");
00798   }
00799 
00807   mfitsio_write_image_execute(fptr, img, 0, info->naxis, lnaxes, fpixels,
00808                               lpixels, &status);
00809 
00810   if (status)
00811   {
00812     MFITSIO_PRINTF ("15 An error occurred: %d\n", status);
00813     MFITSIO_ERR ("Could not write the image to the FITS file.");
00814   }
00815 
00817   fits_close_file (fptr, &status);
00818   if (status)
00819   {
00820     MFITSIO_PRINTF ("39 An error occurred: %d\n", status);
00821     MFITSIO_ERR ("Could not close file.");
00822   }
00823 }
00824 
00825 void mfitsio_write_image_execute(fitsfile *file,
00826                                  mxArray *img, long group, long naxis,
00827                                  long *naxes, long *fpixel, long *lpixel,
00828                                  int *status) {
00829   mxClassID classId;
00830   void *data;
00831   data = mxGetData(img);
00832   classId = mxGetClassID (img);
00833   switch (classId)
00834     {
00835     case mxSINGLE_CLASS:
00836       fits_write_subset_flt(file, group, naxis, naxes, fpixel, lpixel,
00837                             (float*)data, status);
00838       break;
00839     case mxDOUBLE_CLASS:
00841       fits_write_subset_dbl(file, group, naxis, naxes, fpixel, lpixel,
00842                             (double*)data, status);
00843       break;
00844     case mxUINT8_CLASS:
00846       fits_write_subset_byt(file, group, naxis, naxes, fpixel, lpixel,
00847                             (unsigned char*)data, status);
00848       break;
00849     case mxUINT16_CLASS:
00851       fits_write_subset_usht(file, group, naxis, naxes, fpixel, lpixel,
00852                              (unsigned short int*)data, status);
00853       break;
00854     case mxINT16_CLASS:
00856       fits_write_subset_sht(file, group, naxis, naxes, fpixel, lpixel,
00857                             (short int*)data, status);
00858       break;
00859     case mxUINT32_CLASS:
00862       fits_write_subset_uint(file, group, naxis, naxes, fpixel, lpixel,
00863                             (unsigned int*)data, status);
00864       break;
00865     case mxINT32_CLASS:
00867       fits_write_subset_int(file, group, naxis, naxes, fpixel, lpixel,
00868                             (int*)data, status);
00869       break;
00870     case mxLOGICAL_CLASS:
00872       break;
00873     case mxUINT64_CLASS:
00875 #ifdef HAVE_LONGLONG
00876       fits_write_subset_lnglng(file, group, naxis, naxes, fpixel, lpixel,
00877                             (LONGLONG*)data);
00878                             
00879       break;
00880 #endif
00881     case mxINT64_CLASS:
00883 #ifdef HAVE_LONGLONG
00884       fits_write_subset_lnglng(file, group, naxis, naxes, fpixel, lpixel,
00885                             (LONGLONG*)data);
00886       break;
00887 #endif
00888     default:
00890       MFITSIO_PRINTF
00891         ("51 Unsupported data type for image.\n");
00892       MFITSIO_ERR ("Cannot continue.");
00893       break;
00894     }
00895 }
00896 
00904 void
00905 mfitsio_delete_keyword (const char *filename, const char *keyword)
00906 {
00907   fitsfile *fptr = 0;
00908   int status = 0;
00909   mfitsio_record *record = 0;
00910   char *tmp = 0;
00911 
00912   fits_open_file (&fptr, filename, MFITSIO_READWRITE, &status);
00913   if (status)
00914   {
00915     MFITSIO_PRINTF ("16 An error occurred: %d\n", status);
00916     MFITSIO_ERR ("Could not open file for writing.");
00917   }
00918   fits_movabs_hdu (fptr, 1, 0, &status);
00919   if (status)
00920   {
00921     MFITSIO_PRINTF ("17 An error occurred: %d\n", status);
00922     MFITSIO_ERR ("Could not access header.");
00923   }
00924   if (!mfitsio_forbidden (keyword))
00925   {
00926     tmp = strdup (keyword);
00927     fits_delete_key (fptr, tmp, &status);
00928     free (tmp);
00929   }
00930   else
00931   {
00932     MFITSIO_PRINTF ("19 It is forbidden to delete the keyword '%s'\n", keyword);
00933     MFITSIO_ERR ("Could not delete keyword.");
00934   }
00935 
00936   if (status)
00937   {
00938     MFITSIO_PRINTF ("18 An error occurred: %d\n", status);
00939     MFITSIO_ERR ("Could not delete keyword.");
00940   }
00941 }
00942 
00950 void
00951 mfitsio_write_header (const char *filename, const mxArray * header)
00952 {
00953   fitsfile *fptr = 0;
00954   int numFields = 0, i = 0, status = 0, singleton = 0, buflen = 0;
00955   mxClassID classId;
00956   char *fieldName = 0, *stringval = 0;
00957   float floatval = 0.0f;
00958   double dblval = 0.0;
00959   unsigned char ucharval = 0;
00960   char charval = 0;
00961   short shortval = 0;
00962   unsigned short ushortval = 0;
00963   int intval = 0;
00964   unsigned int uintval = 0;
00965   long longval = 0;
00966   unsigned long ulongval = 0;
00967   mxArray *field = 0;
00968 
00969   int naxis = 2;
00970   int bitpix = 8;
00971   long naxes[] = { 1, 1, 1 };
00972 
00975   if (access (filename, F_OK))
00976   {
00977     fits_create_file (&fptr, filename, &status);
00978     if (status)
00979     {
00980       MFITSIO_PRINTF ("20 An error occurred: %d\n", status);
00981       MFITSIO_ERR ("Could not create new FITS file.");
00982     }
00984     fits_create_img (fptr, bitpix, naxis, naxes, &status);
00985     if (status)
00986     {
00987       MFITSIO_PRINTF ("21 An error occurred: %d\n", status);
00988       MFITSIO_ERR ("Could not create an image within the FITS file.");
00989     }
00992     fits_write_key_longwarn (fptr, &status);
00993     if (status)
00994     {
00995       MFITSIO_PRINTF ("22 An error occurred: %d\n", status);
00996       MFITSIO_ERR ("Could not write information to header.");
00997     }
00998   }
00999   else
01000   {
01002     fits_open_file (&fptr, filename, MFITSIO_READWRITE, &status);
01003   }
01004   if (status)
01005   {
01006     MFITSIO_PRINTF ("23 An error occurred: %d\n", status);
01007     MFITSIO_ERR ("Could not open file for writing.");
01008   }
01010   fits_movabs_hdu (fptr, 1, 0, &status);
01011   if (status)
01012   {
01013     MFITSIO_PRINTF ("24 An error occurred: %d\n", status);
01014     MFITSIO_ERR ("Could not access header.");
01015   }
01017   numFields = mxGetNumberOfFields (header);
01018 
01021   for (i = 0; i < numFields; i++)
01022   {
01024     fieldName = (char *) mxGetFieldNameByNumber (header, i);
01025 
01030     if (!mfitsio_forbidden (fieldName))
01031     {
01033       field = mxGetField (header, 0, fieldName);
01034       classId = mxGetClassID (field);
01035 
01038       singleton = mfitsio_is_scalar (field);
01039 
01040       if (singleton || classId == mxCHAR_CLASS)
01041       {
01042 
01044         switch (classId)
01045         {
01046         case mxSINGLE_CLASS:
01048           floatval = (float) mfitsio_get_mscalar (field);
01049           fits_update_key (fptr, TFLOAT, fieldName, &floatval, 0, &status);
01050           break;
01051         case mxDOUBLE_CLASS:
01053           dblval = (double) mfitsio_get_mscalar (field);
01054           fits_update_key (fptr, TDOUBLE, fieldName, &dblval, 0, &status);
01055           break;
01056         case mxCHAR_CLASS:
01058           buflen = (mxGetM (field) * mxGetN (field)) + 1;
01059           stringval = MFITSIO_CALLOC (buflen, sizeof (char));
01060 
01061           status = mxGetString (field, stringval, buflen);
01062           fits_update_key_longstr (fptr, fieldName, stringval, 0, &status);
01063           break;
01064         case mxUINT8_CLASS:
01066           ucharval = (unsigned char) mfitsio_get_mscalar (field);
01067           fits_update_key (fptr, TBYTE, fieldName, &ucharval, 0, &status);
01068           break;
01069         case mxINT8_CLASS:
01071           charval = (char) mfitsio_get_mscalar (field);
01072           fits_update_key (fptr, TBYTE, fieldName, &charval, 0, &status);
01073           break;
01074         case mxUINT16_CLASS:
01077           ushortval = (unsigned short) mfitsio_get_mscalar (field);
01078           fits_update_key (fptr, TUSHORT, fieldName, &ushortval, 0, &status);
01079           break;
01080         case mxINT16_CLASS:
01082           shortval = (short) mfitsio_get_mscalar (field);
01083           fits_update_key (fptr, TSHORT, fieldName, &shortval, 0, &status);
01084           break;
01085         case mxUINT32_CLASS:
01087           ulongval = (unsigned int) mfitsio_get_mscalar (field);
01088           fits_update_key (fptr, TUINT, fieldName, &ulongval, 0, &status);
01089           break;
01090         case mxINT32_CLASS:
01092           longval = (int) mfitsio_get_mscalar (field);
01093           fits_update_key (fptr, TINT, fieldName, &longval, 0, &status);
01094           break;
01095         case mxLOGICAL_CLASS:
01097           intval = (int) mfitsio_get_mscalar (field);
01098           fits_update_key (fptr, TLOGICAL, fieldName, &intval, 0, &status);
01099           break;
01100         default:
01102           MFITSIO_PRINTF
01103             ("25 Unsupported data type for header field \'%s\'.\n", fieldName);
01104           MFITSIO_ERR ("Cannot continue.");
01105           break;
01106         }
01107       }
01108       else
01109       {
01112         MFITSIO_PRINTF
01113           ("26 There exists a non-singleton dimension in the value of \'"
01114            "%s\'. A scalar is required. The field will be ignored.\n");
01115         MFITSIO_WARN ("A header field will not be written.");
01116       }
01118       if (status)
01119       {
01120         MFITSIO_PRINTF
01121           ("27 An error occurred: %d\n. Could not write key for field "
01122            "\'%s\'.", status, fieldName);
01123         MFITSIO_ERR ("Could not write key.");
01124       }
01125     }
01126   }
01128   fits_close_file (fptr, &status);
01129   if (status)
01130   {
01131     MFITSIO_PRINTF ("28 An error occurred: %d\n", status);
01132     MFITSIO_ERR ("Could not create an image within the FITS file.");
01133   }
01134 }
01135 
01144 void
01145 mfitsio_write_info (fitsfile * fptr, const mfitsio_info * info)
01146 {
01147   char naxisf[8];
01148   int status = 0, i = 0, val = 0;
01153   val = info->naxis;
01154   fits_update_key (fptr, TINT, "NAXIS", &val, "", &status);
01155   if (status)
01156   {
01157     MFITSIO_PRINTF ("29 An error occurred: %d\n", status);
01158     MFITSIO_ERR ("Could not store NAXIS information.");
01159   }
01160   val = info->bitpix;
01161   fits_update_key (fptr, TINT, "BITPIX", &val, "", &status);
01162   if (status)
01163   {
01164     MFITSIO_PRINTF ("30 An error occurred: %d\n", status);
01165     MFITSIO_ERR ("Could not store BITPIX information.");
01166   }
01168   for (i = 0; i < info->naxis; i++)
01169   {
01170     snprintf (naxisf, 7, "NAXIS%d", i);
01171     fits_update_key (fptr, TINT, naxisf, (info->naxes + i), "", &status);
01172     if (status)
01173     {
01174       MFITSIO_PRINTF ("31 An error occurred: %d\n", status);
01175       MFITSIO_ERR ("Could not store NAXISX information.");
01176     }
01177   }
01178 }
01179 
01189 mfitsio_info *
01190 mfitsio_calc_info (const mxArray * img)
01191 {
01192   mfitsio_info *retval = 0;
01193   const int *naxesm;
01194   int i;
01195   mxClassID classId;
01197   retval = (mfitsio_info *) MFITSIO_MALLOC (sizeof (mfitsio_info));
01198 
01200   retval->naxis = mxGetNumberOfDimensions (img);
01201   retval->naxes = (int *) MFITSIO_MALLOC (sizeof (int) * retval->naxis);
01202   retval->check = 0;
01203 
01205   naxesm = mxGetDimensions (img);
01206   for (i = 0; i < retval->naxis; i++)
01207   {
01208     retval->naxes[i] = naxesm[i];
01209   }
01210 
01212   classId = mxGetClassID (img);
01213   switch (classId)
01214   {
01215   case mxSINGLE_CLASS:
01216     retval->bitpix = FLOAT_IMG;
01217     break;
01218   case mxDOUBLE_CLASS:
01219     retval->bitpix = DOUBLE_IMG;
01220     break;
01221   case mxUINT8_CLASS:
01222     retval->bitpix = BYTE_IMG;
01223     break;
01224   case mxINT8_CLASS:
01225     retval->bitpix = BYTE_IMG;
01226     break;
01227   case mxUINT16_CLASS:
01228     retval->bitpix = SHORT_IMG;
01229     break;
01230   case mxINT16_CLASS:
01231     retval->bitpix = SHORT_IMG;
01232     break;
01233   case mxUINT32_CLASS:
01234     retval->bitpix = LONG_IMG;
01235     break;
01236   case mxINT32_CLASS:
01237     retval->bitpix = LONG_IMG;
01238     break;
01239   case mxCHAR_CLASS:
01240   case mxUINT64_CLASS:
01241   case mxINT64_CLASS:
01242   case mxLOGICAL_CLASS:
01243   default:
01244     MFITSIO_PRINTF ("32 Unsupported data type for header image.\n");
01245     MFITSIO_ERR ("Cannot continue.");
01246     break;
01247   }
01248   return retval;
01249 }
01250 
01259 int
01260 mfitsio_is_scalar (const mxArray * array)
01261 {
01262   int ndim = 0, nsd = 1, i = 0;
01263   const int *dims;
01265   ndim = mxGetNumberOfDimensions (array);
01266 
01268   dims = mxGetDimensions (array);
01269 
01271   for (i = 0; i < ndim; i++)
01272   {
01273 
01275     if (dims[i] != 1)
01276     {
01277       nsd = 0;
01278       break;
01279     }
01280   }
01281   return nsd;
01282 }
01283 
01296 int
01297 mfitsio_forbidden (const char *s)
01298 {
01299   return !strncmp (s, "NAXIS", 5) || !strcmp (s, "SIMPLE")
01300     || !strcmp (s, "_nil_") || !strcmp (s, "BITPIX");
01301 }
01302 
01315 int
01316 mfitsio_ignore_card (const char *s)
01317 {
01318   return !strncmp (s, "COMMENT  ", 9) || !strncmp (s, "LONGSTRN=", 9)
01319     || !strncmp (s, "EXTEND  =", 9) || !strncmp (s, "SIMPLE  =", 9);
01320 }
01321 
01333 mxArray *
01334 mfitsio_get_mfield (const mxArray * array, int index, const char *field_name)
01335 {
01336   mxArray *field;
01337   if (array == 0)
01338   {
01339     MFITSIO_PRINTF ("34 Header does not exist.\n");
01340     MFITSIO_ERR ("Cannot continue.");
01341   }
01342   field = mxGetField (array, index, field_name);
01343   if (field == 0)
01344   {
01345     MFITSIO_PRINTF ("35 Field '%s' not found in MATLAB struct array.\n",
01346                      field_name);
01347     MFITSIO_ERR ("Cannot continue.");
01348   }
01349   return field;
01350 }
01351 
01363 double
01364 mfitsio_get_mscalar (const mxArray * array)
01365 {
01366   return (array) ? mfitsio_get_mscalar (array) : 0;
01367 }
01368 
01380 void mfitsio_check_coordinate(const mxArray *crd1, const mxArray *crd2,
01381                               const int naxis, const int *naxes) {
01382   int i;
01383   double *data1 = 0, *data2 = 0;
01384   if (mxGetNumberOfDimensions(crd1) > 2) {
01385     MFITSIO_PRINTF ("40 Too many  in your starting pixel coordinate."
01386                     " It must be a 1x%d vector.\n", naxis);
01387     MFITSIO_ERR ("Invalid coordinate.");
01388   }
01389   if (mxGetNumberOfDimensions(crd2) > 2) {
01390     MFITSIO_PRINTF ("41 Too many dimensions in your ending pixel coordinate. "
01391                     "It must be a 1x%d vector.\n", naxis);
01392     MFITSIO_ERR ("Invalid coordinate.");
01393   }
01394   if (mxGetM(crd1) != 1) {
01395     MFITSIO_PRINTF ("42 Too many rows in your starting pixel coordinate. "
01396                     "It must be a 1x%d vector.\n", naxis);
01397     MFITSIO_ERR ("Invalid coordinate.");
01398   }
01399   if (mxGetM(crd2) != 1) {
01400     MFITSIO_PRINTF ("43 Too many rows in your ending pixel coordinate. "
01401                     "It must be a 1x%d vector.\n", naxis);
01402     MFITSIO_ERR ("Invalid coordinate.");
01403   }
01404   if (mxGetN(crd1)<naxis) {
01405     MFITSIO_PRINTF ("44 Not enough values in row vector defining"
01406                     " starting pixel coordinate. "
01407                     "It must be a 1x%d vector.\n", naxis);
01408     MFITSIO_ERR ("Invalid coordinate.");
01409   }
01410   if (mxGetN(crd1)>naxis) {
01411     MFITSIO_PRINTF ("45 Too many values in row vector defining"
01412                     " starting pixel coordinate. "
01413                     "It must be a 1x%d vector.\n", naxis);
01414     MFITSIO_ERR ("Invalid coordinate.");
01415   }
01416   if (mxGetN(crd2)<naxis) {
01417     MFITSIO_PRINTF ("46 Not enough values in row vector defining"
01418                     " end pixel coordinate."
01419                     "It must be a 1x%d vector.\n", naxis);
01420     MFITSIO_ERR ("Invalid coordinate.");
01421   }
01422   if (mxGetN(crd2)>naxis) {
01423     MFITSIO_PRINTF ("47 Too many values in row vector defining"
01424                     " starting pixel coordinate."
01425                     "It must be a 1x%d vector.\n", naxis);
01426     MFITSIO_ERR ("Invalid coordinate.");
01427   }
01428   data1 = mxGetData(crd1);
01429   data2 = mxGetData(crd2);
01430   for (i = 0; i < naxis; i++ ) {
01431     if (data1[ i ] < 1 || data1[ i ] > naxes[ i ]) {
01432       MFITSIO_PRINTF("36 Invalid component at index %d of starting"
01433                      " pixel: %d\n", i + 1, (int)data1[ i ]);
01434       MFITSIO_ERR("Starting coordinate out of range.");
01435     }
01436     if (data2[ i ] < 1 || data2[ i ] > naxes[ i ]) {
01437       MFITSIO_PRINTF("37 Invalid component at index %d of ending"
01438                      " pixel: %d\n", i + 1, (int)data2[ i ]);
01439       MFITSIO_ERR("Ending coordinate out of range.");
01440     }
01441     if (data2[ i ] < data1[ i ]) {
01442       MFITSIO_PRINTF("38 Component at index %d of starting pixel must be
01443                       less than or equal to component %d of ending pixel.\n");
01444       MFITSIO_ERR("Invalid subset (hyper)region.");
01445     }
01446   }
01447 }
01448 
01456 void mfitsio_check_size_info(const mxArray *size_unit) {
01457 
01458   int i;
01459   double *data = 0;
01460   int naxis;
01461   if (mxGetNumberOfDimensions(size_unit) > 2) {
01462     MFITSIO_PRINTF ("48 Too many dimensions defining your axes."
01463                     " The size unit must be a 1xN vector.\n");
01464     MFITSIO_ERR ("Invalid coordinate.");
01465   }
01466 
01467   if (mxGetM(size_unit) != 1) {
01468     MFITSIO_PRINTF ("49 Too many rows defining your axes. "
01469                     " The size unit must be a 1xN vector.\n", naxis);
01470     MFITSIO_ERR ("Invalid coordinate.");
01471   }
01472   data = mxGetData(size_unit);
01473   naxis = mxGetN(size_unit);
01474   for (i = 0; i < naxis; i++ ) {
01475     if (data[ i ] < 0) {
01476       MFITSIO_PRINTF("50 Invalid component at index %d of size unit"
01477                      ": %d\n", i + 1, (int)data[ i ]);
01478       MFITSIO_ERR("Size unit out of range.");
01479     }
01480   }
01481 }
01482 
01483 
01493 long *mfitsio_convert_dbl2long(const double *dbl, const int size) {
01494   int i;
01495   long *l = (long*)MFITSIO_MALLOC(sizeof(long)*size);
01496   for (i = 0; i < size; i++ ) {
01497     l[ i ] = (long)dbl[ i ];
01498   }
01499   return l;
01500 }
01501 
01511 long *mfitsio_convert_int2long(const int *in, const int size) {
01512   int i;
01513   long *l = (long*)MFITSIO_MALLOC(sizeof(long)*size);
01514   for (i = 0; i < size; i++ ) {
01515     l[ i ] = (long)in[ i ];
01516   }
01517   return l;
01518 }
01519 
01528 long *mfitsio_create_ones_vector(const int size) {
01529   int i;
01530   long *l = (long*)MFITSIO_MALLOC(sizeof(long)*size);
01531   for (i = 0; i < size; i++ ) {
01532     l[ i ] = 1;
01533   }
01534   return l;
01535 }
01536 
01547 int *mfitsio_get_region_size(const long *crd1, const long *crd2,
01548                              const int naxis) {
01549   int i;
01550   int *diff = 0;
01551   diff = (int*)MFITSIO_MALLOC(sizeof(int)*naxis);
01552   for (i = 0; i < naxis; i++ ) {
01553     diff[ i ] = (int)(crd2[ i ] - crd1[ i ]) + 1;
01554   }
01555   return diff;
01556 }
01557 
01568 mxArray *mfitsio_get_lpixels(const mxArray *spixels,
01569                              const mxArray *img, 
01570                              int naxis) {
01571   int i;
01572   int dim[] = { 1, 1 };
01573   mxArray *array;
01574   double *sdata;
01575   double *ldata;
01576   const int *imgdims;
01577   int nimgdims;
01578   dim[ 1 ] = naxis;
01579   imgdims = mxGetDimensions(img);
01580   array = mxCreateNumericArray (2, dim, mxDOUBLE_CLASS, mxREAL);  
01581   sdata = mxGetData(spixels);
01582   ldata = mxGetData(array);
01583   nimgdims = mxGetNumberOfDimensions(img);
01584   for (i = 0; i < nimgdims; i++ ) {
01585     ldata[ i ] = sdata[ i ] + (double)imgdims[ i ] - 1;
01586   }
01587   for (; i < naxis; i++ ) {
01588     ldata[ i ] = sdata[ i ];
01589   }
01590   return array;
01591 }
01592 
01593 char *mfitsio_logical_m2f(const mxArray *img) {
01594   const int *imgdims;
01595   mxLogical *logicals;
01596   char *retval = 0;
01597   long prod = 1, j;
01598   int nimgdims, i;
01599   logicals = mxGetLogicals(img);
01600   imgdims = mxGetDimensions(img);
01601   nimgdims = mxGetNumberOfDimensions(img);
01602   for (i = 0; i < nimgdims; i++ ) {
01603     prod *= ((long)imgdims[ i ]);
01604   }
01605   retval = (char*)MFITSIO_MALLOC(sizeof(char) * prod);
01606   for (j = 0; j < prod; j++ ) {
01607     retval[ j ] = logicals[ j ] ? 'T' : 'F';
01608   }
01609   return retval;
01610 }