my_eclipse.c

00001 /* 
00002 #include <stdio.h>
00003 #include <stdlib.h>
00004 #include <string.h>
00005 #include <unistd.h>
00006 */
00007 #include <stdarg.h>
00008 #include "my_eclipse.h"
00009 #include "my_matrix.h"
00010 #include "my_ipow.h"
00011 
00012 #define ALLOC_CELL_NOALLOC      -1
00013 #define ALLOC_CELL_RAM           0
00014 #define ALLOC_CELL_SWAP          1
00015 
00016 #define SRCFILENAME_SZ          80
00017 #define TMPFILENAME_SZ          256
00018 
00019 #define CLEQSIZE              65534 /* limit class equivalence table size */
00020 
00021 #define LOGNAMESZ                       512
00022 #define HOSTNAMESZ                      512
00023 #define TIMESTRSZ                       512
00024 #define  ARC_RANGE_FACT          3.0
00025 #define  TRESH_MEDIAN_MIN        0.0
00026 #define  TRESH_SIGMA_MAX         200.0
00027 #define dtiny(a) ((a)<0?(a)> -1.e-30:(a)<1.e-30)
00028 
00029 static  int             file_reg = 0 ;
00030 
00031 static int select_valid_arcs(int, int, cc_stats_t *, int, int *, int **,
00032                              int **, orientation_t) ;
00033 
00034 /*
00035  * The 'keytuple' type is strictly internal to this module.
00036  * It describes FITS cards as tuples (key,value,comment,line), where key
00037  * is always a non-NULL character string, value and comment are
00038  * allowed to be NULL. 'line' is a string containing the line as it
00039  * has been read from the input FITS file (raw). It is set to NULL if
00040  * the card is modified later. This allows in output two options:
00041  * either reconstruct the FITS lines by printing key = value / comment
00042  * in a FITS-compliant way, or output the lines as they were found in
00043  * input, except for the modified ones.
00044  *
00045  * The following functions are associated methods
00046  * to this data structure:
00047  * keytuple_new()               constructor
00048  * keytuple_del()               destructor
00049  * keytuple_cmp()               comparison function based on key only
00050  * keytuple_dmp()               dumps a keytuple to stdout
00051  */
00052 
00053 
00054 typedef struct _keytuple_ {
00055         char    *       key ;
00056         char    *       val ;
00057         char    *       com ;
00058         char    *       lin ;
00059 } keytuple ;
00060 
00061 
00062 
00063 
00064 
00065 static keytuple *       keytuple_new(char *, char *, char *, char *);
00066 static void             keytuple_del(keytuple * k);
00067 static int                      keytuple_cmp(const void * key1, const void * key2);
00068 static void             keytuple_dmp(keytuple * k);
00069 
00070 
00071 
00072 
00073 
00074 
00075 static int is_blank_line(char * s)
00076 {
00077         int             i ;
00078 
00079         for (i=0 ; i<strlen(s) ; i++) {
00080                 if (s[i]!=' ') return 0 ;
00081         }
00082         return 1 ;
00083 }
00084 
00085 static int gauss_pivot(double *ptra, double *ptrc, int n);
00086 
00087 
00088 /*
00089  * The following parameters set up the limits for this memory module.
00090  * They are defined as global variables to this module and should not be
00091  * seen outside except through member functions get/set_memory_parameter().
00092  */
00093 
00094 static size_t emem_ramlimit                     = 1024 ;        /* megabytes */
00095 static size_t emem_vm_limit                     = 1024 ;        /* megabytes */
00096 static size_t emem_vm_pagesize          = 1024 ;        /* kilobytes */
00097 static char   emem_tmpdirname[1024]     = "." ;
00098 
00099 
00100 
00101 static void * ram_malloc(size_t, char *, int);
00102 static void * swap_malloc_x(size_t, char *, int);
00103 /*
00104 static void * ram_calloc(size_t, size_t, char *, int);
00105 static void   ram_free(void *, char *, int);
00106 static void   swap_free_x(void *, char *, int);
00107 static void   cleanup_mess(void);
00108 static void   cleanup_mess_wrapper(void);
00109 static int    signal_catch(int sig, void (*f)());
00110 static void   fatal_signal_handler(void);
00111 */
00112 
00113 
00114 /* This structure holds all allocated pointers (forward linked list)    */
00115 typedef struct _MALLOC_CELL_ {
00116         void                            *       pointer;
00117         size_t                                  size;
00118         char                                    filename[SRCFILENAME_SZ] ;
00119         int                                             lineno ;
00120         int                                             flavour ;
00121         char                                    swapfilename[TMPFILENAME_SZ] ;
00122         int                                             swapfd ;
00123         struct _MALLOC_CELL_*   next ;
00124         struct _MALLOC_CELL_*   prev ;
00125 } malloc_cell;
00126 
00127 
00128 /* This structure holds information about everything allocated so far   */
00129 static malloc_cell * malloc_table = NULL ;
00130 
00131 /*
00132  * Global variables. Their visibility is limited to this module.
00133  * Access and modification of their values is done through associated
00134  * member function calls.
00135  */
00136 static  size_t  total_allocated_RAM = 0 ;
00137 static  size_t  total_allocated_SWAP = 0 ;
00138 static  int             nswapfiles = 0 ;
00139 static  int             nalloc_ptrs = 0 ;
00140 
00141 
00142 
00143 int eclipse_verbose = 0 ;
00144 int eclipse_debug   = 0 ;
00145 int eclipse_logfile = 0 ;
00146 
00147 static char eclipse_logfilename[LOGNAMESZ];
00148 /* Chop a string: private function */
00149 static void chop_r(char * s)
00150 {
00151         if (s) {
00152                 if (strlen(s)>1) {
00153                         while (s[strlen(s)-1] == (char)'\n') s[strlen(s)-1]=(char)0 ;
00154                 }
00155         }
00156 }
00157 
00158 
00159 
00160 /*-------------------------------------------------------------------------*/
00174 /*--------------------------------------------------------------------------*/
00175 
00176 OneImage *
00177 new_image(
00178     int     size_x,
00179     int     size_y
00180 )
00181 {
00182     OneImage    *       new_image ;
00183 
00184     if ((size_x<1) ||
00185         (size_x>MAX_COLUMN_NUMBER) ||
00186         (size_y<1) ||
00187         (size_y>MAX_LINE_NUMBER)) {
00188         e_error("requested size: [%dx%d]", size_x, size_y) ;
00189         e_error("cannot create such size: aborting image creation") ;
00190         return NULL ;
00191     }
00192 
00193     new_image = calloc(1, sizeof(OneImage)) ;
00194 
00195     new_image->lx = size_x ;
00196     new_image->ly = size_y ;
00197     new_image->nbpix = size_x * size_y ;
00198     new_image->data = calloc(size_x * size_y, sizeof(pixelvalue));
00199         new_image->reg.active = 0 ;
00200     return new_image ;
00201 }
00202 
00203 /*-------------------------------------------------------------------------*/
00213 /*--------------------------------------------------------------------------*/
00214 
00215 OneImage *
00216 load_image(
00217     const char *filename
00218 )
00219 {
00220     OneImage    *loaded ;
00221     cube_info    *fileinfo ;
00222 
00223     fileinfo = get_cube_info(filename) ;
00224     if (fileinfo == NULL) {
00225         e_error("in file [%s]: aborting cube loading", filename) ;
00226         return NULL ;
00227     }
00228 
00229     /* Load first plane in cube only    */
00230     loaded = load_plane_from_fits(filename, fileinfo, 0) ;
00231     free(fileinfo) ;
00232     if (loaded == NULL) {
00233         e_error("cannot load file [%s]", filename) ;
00234     }
00235 
00236     return loaded ;
00237 }
00238 
00239 
00240 
00241 
00242 
00243 /*-------------------------------------------------------------------------*/
00255 /*--------------------------------------------------------------------------*/
00256 
00257 cube_info *
00258 get_cube_info(char *filename)
00259 {
00260         char            *       value ;
00261     int                 naxes ;
00262     cube_info    *      fileinfo ;
00263 
00264     /* Check file existence */
00265     if (file_exists(filename) != 1) {
00266         e_error("file %s not found", filename) ;
00267         return((cube_info*)NULL) ;
00268     }
00269 
00270     /* Allocate space for returned info */
00271     fileinfo = malloc(sizeof(cube_info)) ;
00272 
00273     /* Get the header size  */
00274     fileinfo->headersize = fits_get_header_size(filename);
00275         if (fileinfo->headersize<1) {
00276                 free(fileinfo);
00277                 return NULL ;
00278         }
00279 
00280     /* find the number of axes  */
00281 
00282         naxes = 0 ;
00283         value = query_fits_header(filename, "NAXIS");
00284         if (value!=NULL) {
00285                 sscanf(value, "%d", &naxes);
00286         } else {
00287                 e_error("missing key in header: NAXIS");
00288                 free(fileinfo);
00289                 return NULL ;
00290         }
00291 
00292     /*
00293      * Cubes may not have less than one dimension or more
00294      * than 3 dimensions.
00295      * One dimensional cubes are processed as single-lined images.
00296      */
00297 
00298     /* Get number of axes */
00299     if  ((naxes < 1) || (naxes > 3))    {
00300         e_error("cannot handle cube with %d axes", naxes) ;
00301                 free(fileinfo);
00302         return NULL ;
00303     }
00304 
00305     /* Then read the dimension value for each axis  */
00306     /* NAXIS1 is always present */
00307         fileinfo->lx = 0 ;
00308         value = query_fits_header(filename, "NAXIS1");
00309         if (value!=NULL) {
00310                 sscanf(value, "%d", &(fileinfo->lx)) ;
00311         } else {
00312                 e_error("missing key in header: NAXIS1");
00313                 free(fileinfo);
00314         return NULL ;
00315         }
00316 
00317     if ((fileinfo->lx < 1) || (fileinfo->lx > MAX_COLUMN_NUMBER)) {
00318         e_error("cannot process cube with NAXIS1 = %d", fileinfo->lx) ;
00319                 free(fileinfo);
00320         return NULL ;
00321     }
00322 
00323     /*
00324      * Now get second axis size value.
00325      */
00326 
00327     if (naxes == 1) {
00328         fileinfo->ly = 1 ;
00329     } else {
00330         /* Read value and assign ly */
00331         fileinfo->ly = 0 ;
00332                 value = query_fits_header(filename, "NAXIS2");
00333                 if (value!=NULL) {
00334                         sscanf(value, "%d", &(fileinfo->ly)) ;
00335                 } else {
00336                         e_error("missing key in header: NAXIS2");
00337                         free(fileinfo);
00338                         return NULL ;
00339                 }
00340     }
00341 
00342     if ((fileinfo->ly < 1) || (fileinfo->ly > MAX_LINE_NUMBER)) {
00343         e_error("cannot process cube with NAXIS2 = %d", fileinfo->ly) ;
00344                 free(fileinfo);
00345         return NULL ;
00346     }
00347 
00348     fileinfo->timeflag = ((fileinfo->ly == 129) || (fileinfo->ly == 257)) ?
00349                             1 :
00350                             0 ;
00351 
00352 
00353     /* For 2D cubes, this is all we have to do      */
00354     /* For 3D cubes, get the third dimension size   */
00355    if (naxes <= 2) {
00356         fileinfo->n_im = 1 ;
00357     } else /* It is a 3D cube   */ {
00358                 fileinfo->n_im = 0 ;
00359                 value = query_fits_header(filename, "NAXIS3");
00360                 if (value!=NULL) {
00361                         sscanf(value, "%d", &(fileinfo->n_im)) ;
00362                 } else {
00363                         e_error("missing key in header: NAXIS3");
00364                         free(fileinfo);
00365                         return NULL ;
00366                 }
00367 
00368         if ((fileinfo->n_im < 1) || (fileinfo->n_im > MAX_IMAGE_NUMBER)) {
00369             e_error("cannot process cube with NAXIS3 = %d", fileinfo->n_im) ;
00370             free(fileinfo) ;
00371             return NULL ;
00372         }
00373         /* Third axis size : OK */
00374     }
00375 
00376     /* Test data type   */
00377 
00378     fileinfo->ptype = 0 ;
00379         value = query_fits_header(filename, "BITPIX");
00380         if (value!=NULL) {
00381                 sscanf(value, "%d", &(fileinfo->ptype)) ;
00382         } else {
00383                 e_error("missing key in header: BITPIX");
00384                 free(fileinfo);
00385                 return NULL ;
00386         }
00387         if (BYTESPERPIXEL(fileinfo->ptype) == 0) {
00388         e_error("cannot process cube with resolution of %s", value) ;
00389         free(fileinfo) ;
00390         return NULL ;
00391     }
00392 
00393     /* Get BSCALE keyword, or assign 1.0 if not present */
00394 
00395         fileinfo->b_scale = 1.0 ;
00396         value = query_fits_header(filename, "BSCALE");
00397         if (value!=NULL) {
00398                 sscanf(value, "%lg", &(fileinfo->b_scale)) ;
00399                 /* issue a warning is BSCALE=0 */
00400                 if (atoi(value)==0) {
00401                         e_warning("BSCALE=0 in header: using 1");
00402                         fileinfo->b_scale = 1.0 ;
00403                 }
00404         }
00405 
00406     /* Get BZERO keyword, or assign 0.0 if not present  */
00407 
00408         fileinfo->b_zero = 0.0 ;
00409         value = query_fits_header(filename, "BZERO");
00410         if (value!=NULL)
00411                 sscanf(value, "%lg", &(fileinfo->b_zero)) ;
00412 
00413     /* More tests should be added here to test the real validity of     */
00414     /* an Adonis FITS file.                                             */
00415     /* Most of all, no keyword is tested.                               */
00416     /* We only check what is necessary to load data into memory.        */
00417 
00418     return fileinfo ;
00419 }
00420 
00421 
00422 /*-------------------------------------------------------------------------*/
00433 /*--------------------------------------------------------------------------*/
00434 
00435 int file_exists(char * filename)
00436 {
00437         struct stat buf ;
00438         int                     exists ;
00439 
00440         if (stat(filename, &buf)==-1) {
00441                 exists=0;
00442         } else {
00443                 exists=1;
00444         }
00445         return exists;
00446 }
00447 
00448 
00449 /*-------------------------------------------------------------------------*/
00460 /*--------------------------------------------------------------------------*/
00461 
00462 
00463 #define FITS_BLSZ 2880
00464 size_t fits_get_header_size(char * name)
00465 {
00466     FILE    * in ;
00467     char      line[81] ;
00468     int       found = 0 ;
00469     int       count ;
00470     size_t    hs ;
00471 
00472     if ((in = fopen(name, "r")) == NULL) {
00473         fprintf(stderr, "cannot open %s: aborting\n", name) ;
00474         return 0 ;
00475     }
00476     count = 0 ;
00477     while (!found) {
00478                 if (fread(line, 1, 80, in)!=80) {
00479                         break ;
00480                 }
00481         count ++ ;
00482         if (!strncmp(line, "END ", 4)) {
00483             found = 1 ;
00484         }
00485     }
00486         fclose(in);
00487 
00488         if (!found) return 0 ;
00489     /*
00490      * The size is the number of found cards times 80,
00491      * rounded to the closest higher multiple of 2880.
00492      */
00493     hs = (size_t)(count * 80) ;
00494     if ((hs % FITS_BLSZ) != 0) {
00495         hs = (1+(hs/FITS_BLSZ)) * FITS_BLSZ ;
00496     }
00497     return hs ;
00498 }
00499 #undef FITS_BLSZ
00500 
00501 
00502 /*-------------------------------------------------------------------------*/
00527 /*--------------------------------------------------------------------------*/
00528 
00529 
00530 char *
00531 query_fits_header(char * filename, char * keyword)
00532 {
00533     mmap_file   *   mm ;
00534     char        *   where ;
00535     char        *   start ;
00536     char        *   value ;
00537     char            test1, test2 ;
00538         int                             i ;
00539         int                             len ;
00540         int                             different ;
00541         struct stat             fileinfo ;
00542         long int                bufcount;
00543 
00544     mm = mmap_open(filename) ;
00545     if (mm==NULL) return NULL ;
00546         stat(filename, &fileinfo);
00547         bufcount=0;
00548     start = mm->buf ;
00549         where = start ;
00550         len   = (int)strlen(keyword) ;
00551     while (1) {
00552                 different=0 ;
00553                 for (i=0 ; i<len ; i++) {
00554                         if (where[i]!=keyword[i]) {
00555                                 different++ ;
00556                                 break ;
00557                         }
00558                 }
00559                 if (!different) {
00560                         /* Get the 2 characters after the keyword */
00561                         test1 = where[len] ;
00562                         test2 = where[len+1] ;
00563                         /* If first subsequent character is the equal sign, bingo. */
00564                         if (test1=='=') break ;
00565                         /*
00566                          * If first subsequent character is a blank and the one after is
00567                          * equal or another blank, bingo.
00568                          */
00569                         if (test1==' ' && (test2=='=' || test2==' ')) break ;
00570                 }
00571                 if ((where[0]=='E') &&
00572                         (where[1]=='N') &&
00573                         (where[2]=='D') &&
00574                         (where[3]==' ')) {
00575                         /* End of header detected */
00576                         mmap_close(mm);
00577                         return NULL ;
00578                 }
00579                 /* Go to next line in header */
00580                 where += 80 ;
00581                 bufcount +=80;
00582                 if (bufcount>=fileinfo.st_size){
00583                         /* file damaged or not a fits file, avoid segementation
00584                          * fault */
00585                         mmap_close(mm);
00586                         return NULL ;
00587                 }
00588 
00589     }
00590     /* Found the keyword, now get its value */
00591     value = fits_getvalue(where);
00592     mmap_close(mm);
00593     return value;
00594 }
00595 /*-------------------------------------------------------------------------*/
00617 /*--------------------------------------------------------------------------*/
00618 
00619 mmap_file *
00620 mmap_open(char * filename)
00621 {
00622         mmap_file       *       mm ;
00623         char            *       buf ;
00624         int                             fd ;
00625         struct stat             fileinfo ;
00626 
00627         /* Find out file size in bytes (POSIX compliant) */
00628 
00629         if (stat(filename, &fileinfo)!=0) {
00630                 fprintf(stderr, "error getting stat for file %s\n", filename);
00631                 return NULL ;
00632     }
00633         if (fileinfo.st_size<1) return NULL ;
00634 
00635         if ((fd = open(filename, O_RDONLY)) == -1) {
00636                 fprintf(stderr,"cannot open file %s: mmap failed\n", filename) ;
00637                 return NULL ;
00638         }
00639         buf = (char*)mmap(0,
00640                                           fileinfo.st_size,
00641                                           PROT_READ | PROT_WRITE,
00642                                           MAP_PRIVATE,
00643                                           fd,
00644                                           0) ;
00645 
00646         if (buf==(char*)-1) { perror("mmap") ; close(fd) ; return NULL ; }
00647         close(fd) ;
00648 
00649         mm = malloc(sizeof(mmap_file)) ;
00650         mm->fd   = fd ;
00651         mm->buf  = buf ;
00652         mm->size = fileinfo.st_size ;
00653        (void)strcpy(mm->filename, filename);
00654         return mm ;
00655 }
00656 
00657 /*-------------------------------------------------------------------------*/
00670 /*--------------------------------------------------------------------------*/
00671 
00672 
00673 char * fits_getvalue(char * line)
00674 {
00675     static char value[81] ;
00676     int     i ;
00677     int     from, to ;
00678     int     inq ;
00679 
00680         if (line==NULL) {
00681 #ifdef DEBUG_FITSHEADER
00682                 printf("fits_getvalue: NULL input line\n");
00683 #endif
00684                 return NULL ;
00685         }
00686 
00687         /* Special cases */
00688 
00689         /* END has no associated value */
00690         if (!strncmp(line, "END ", 4)) {
00691                 return NULL ;
00692         }
00693         /*
00694          * HISTORY has for value everything else on the line, stripping
00695          * blanks before and after. Blank HISTORY is also accepted.
00696          */
00697         memset(value, 0, 81);
00698 
00699         if (!strncmp(line, "HISTORY ", 8) || !strncmp(line, "        ", 8)) {
00700                 i=7 ;
00701                 /* Strip blanks from the left side */
00702                 while (line[i]==' ' && i<80) i++ ;
00703                 if (i>=80) return NULL ; /* Blank HISTORY */
00704                 from=i ;
00705 
00706                 /* Strip blanks from the right side */
00707                 to=79 ;
00708                 while (line[to]==' ') to-- ;
00709                 /* Copy relevant characters into output buffer */
00710                 strncpy(value, line+from, to-from+1);
00711                 /* Null-terminate the string */
00712                 value[to-from+1] = (char)0;
00713                 return value ;
00714         } else if (!strncmp(line, "COMMENT ", 8)) {
00715         /*
00716          * COMMENT is like HISTORY
00717          */
00718                 /* Strip blanks from the left side */
00719                 i=7 ;
00720                 while (line[i]==' ' && i<80) i++ ;
00721                 if (i>=80) return NULL ;
00722                 from=i ;
00723 
00724                 /* Strip blanks from the right side */
00725                 to=79 ;
00726                 while (line[to]==' ') to-- ;
00727 
00728                 if (to<from) {
00729 #ifdef DEBUG_FITSHEADER
00730                         printf("fits_getvalue: inconsistent value search in COMMENT\n");
00731 #endif
00732                         return NULL ;
00733                 }
00734                 /* Copy relevant characters into output buffer */
00735                 strncpy(value, line+from, to-from+1);
00736                 /* Null-terminate the string */
00737                 value[to-from+1] = (char)0;
00738                 return value ;
00739         }
00740         /*
00741          * General case
00742          * Get past the keyword
00743          */
00744     i=0 ;
00745     while (line[i]!='=' && i<80) i++ ;
00746     if (i>80) {
00747 #ifdef DEBUG_FITSHEADER
00748                 printf("fits_getvalue: no equal sign found on line\n");
00749 #endif
00750                 return NULL ;
00751         }
00752     i++ ;
00753     while (line[i]==' ' && i<80) i++ ;
00754     if (i>80) {
00755 #ifdef DEBUG_FITSHEADER
00756                 printf("fits_getvalue: no value past the equal sign\n");
00757 #endif
00758                 return NULL ;
00759         }
00760     from=i;
00761 
00762         /*
00763          * Now in the value section
00764          * Look for the first slash '/' outside of a string
00765          */
00766     inq = 0 ;
00767     while (i<80) {
00768         if (line[i]=='\'')
00769                         inq=!inq ;
00770         if (line[i]=='/')
00771            if (!inq)
00772                 break ;
00773         i++;
00774     }
00775     i-- ;
00776 
00777         /* Backtrack on blanks */
00778     while (line[i]==' ' && i>=0) i-- ;
00779     if (i<0) {
00780 #ifdef DEBUG_FITSHEADER
00781                 printf("fits_getvalue: error backtracking on blanks\n");
00782 #endif
00783                 return NULL ;
00784         }
00785     to=i ;
00786 
00787         if (to<from) {
00788 #ifdef DEBUG_FITSHEADER
00789                 printf("fits_getvalue: from>to?\n");
00790                 printf("line=[%s]\n", line);
00791 #endif
00792                 return NULL ;
00793         }
00794         /* Copy relevant characters into output buffer */
00795     strncpy(value, line+from, to-from+1);
00796         /* Null-terminate the string */
00797     value[to-from+1] = (char)0;
00798         /*
00799          * Make it pretty: remove head and tail quote, change double
00800          * quotes to simple ones.
00801          */
00802     if (value[0]=='\'')
00803         strcpy(value, fits_pretty_string(value));
00804     return value ;
00805 }
00806 
00807 /*-------------------------------------------------------------------------*/
00817 /*--------------------------------------------------------------------------*/
00818 
00819 
00820 void mmap_close(mmap_file * mm)
00821 {
00822         if (mm==NULL) return ;
00823         if (munmap(mm->buf, mm->size)!=0) perror("munmap") ;
00824         free(mm);
00825         return ;
00826 }
00827 
00828 /*-------------------------------------------------------------------------*/
00843 /*--------------------------------------------------------------------------*/
00844 
00845 OneImage *
00846 load_plane_from_fits(
00847     char        *       filename,
00848     cube_info   *       fileinfo,
00849     int                 n_plane
00850 )
00851 {
00852     OneImage    *       loaded_plane ;
00853     FILE        *       input_file ;
00854     size_t              offset ;
00855     int                 bytes_per_pixel ;
00856     BYTE        *       byte_buffer ;
00857     int                 imsize;
00858 
00859     /* Error handling : check entries   */
00860     if (filename == NULL) return NULL ;
00861     if (fileinfo == NULL) return NULL ;
00862     if ((n_plane<0) || (n_plane>=fileinfo->n_im)) {
00863         e_error("out of bounds: requested plane %d, max is %d",
00864                                 n_plane, fileinfo->n_im) ;
00865         return NULL ;
00866     }
00867 
00868     bytes_per_pixel = BYTESPERPIXEL(fileinfo->ptype) ;
00869 
00870     /* Allocate one plane in bytes  */
00871     byte_buffer = malloc(fileinfo->lx *
00872                          fileinfo->ly *
00873                          bytes_per_pixel * sizeof(BYTE)) ;
00874 
00875     if ((input_file = fopen(filename, "r")) == NULL) {
00876         e_error("Problems in reading %s : aborting extraction", filename);
00877         free(byte_buffer) ;
00878         return NULL ;
00879     }
00880     /*
00881      * The official formula for the offset in file is :
00882      * offset =   header_size +
00883          *                        size_x * size_y * bytes_per_pixel * requested_plane
00884      */
00885 
00886     offset =    (size_t)n_plane *
00887                 (size_t)fileinfo->lx *
00888                 (size_t)fileinfo->ly *
00889                 (size_t)bytes_per_pixel + fileinfo->headersize ;
00890 
00891     fseek(input_file, offset, SEEK_SET) ;
00892     imsize =    fileinfo->lx * fileinfo->ly * bytes_per_pixel ;
00893     if ((fread( byte_buffer,
00894                                 sizeof(BYTE),
00895                                 imsize,
00896                                 input_file))!=imsize) {
00897         e_error("while reading %ld bytes in file %s : offset %ld",
00898                 imsize,
00899                                 filename,
00900                                 offset) ;
00901         free(byte_buffer) ;
00902         fclose(input_file) ;
00903         return NULL ;
00904     }
00905     fclose(input_file) ;
00906 
00907     /* Now create OneImage to receive data in internal format    */
00908     loaded_plane = new_image(fileinfo->lx, fileinfo->ly) ;
00909     if (loaded_plane == NULL) {
00910         e_error("cannot allocate output image: aborting extraction") ;
00911         free(byte_buffer) ;
00912         return NULL ;
00913     }
00914 
00915     convert_fits_to_local(loaded_plane->data,
00916                           byte_buffer,
00917                           loaded_plane->nbpix,
00918                           fileinfo->ptype,
00919                           fileinfo->b_scale,
00920                           fileinfo->b_zero) ;
00921     free(byte_buffer) ;
00922     return loaded_plane ;
00923 }
00924 
00925 
00926 
00927 /*
00928  * ---------------------------- IMPORTANT -----------------------------------
00929  *
00930  *
00931  *      Portability is an issue in this module !
00932  *      The FITS standard assumes the following:
00933  *
00934  *      type            # of bits
00935  *
00936  *      char            8
00937  *      short           16
00938  *      long            32
00939  *      float           32
00940  *      double          64
00941  *
00942  *      All these types are converted to the local internal float representation
00943  *      thus bringing portability issues. To handle this:
00944  *
00945  *      - FITS types are assumed, i.e. the above sizes are hardcoded,
00946  *        instead of using sizeof, the numerical values are used (8,16,32,64)
00947  *      - Only Intel and Motorola endian schemes shall be supported in this
00948  *    version. Other schemes require much smarter algorithms and prevent
00949  *    any optimization in these loops, increasing I/O times a lot.
00950  *  - All types have been redefined, indicating bit numbers. It is
00951  *        therefore COMPULSORY to use these types and no char, short, long
00952  *    whenever the number of bits may be an issue.
00953  */
00954 
00955 /*
00956  * The INTEL/MOTOROLA conversion functions are only declared and compiled
00957  * if the local machine byte ordering scheme is INTEL.
00958  * Even then, they are static (private) to this module.
00959  */
00960 
00961 
00962 #if (BYTE_ORDERING_SCHEME == INTEL)
00963 
00964 /*-------------------------------------------------------------------------*/
00974 /*--------------------------------------------------------------------------*/
00975 
00976 static short16
00977 swap_2_bytes(short16 w)
00978 {
00979         register short16 tmp ;
00980         tmp = (w & (ushort16)0x00FF) ;
00981         tmp = ((w & (ushort16)0xFF00) >> (ushort16)0x08) | (tmp << (ushort16)0x08);
00982         return(tmp) ;
00983 }
00984 
00985 #endif
00986 
00987 
00988 
00989 
00990 
00991 #if (BYTE_ORDERING_SCHEME == INTEL)
00992 
00993 
00994 /*-------------------------------------------------------------------------*/
01006 /*--------------------------------------------------------------------------*/
01007 
01008 static long32
01009 swap_4_bytes(long32 dw)
01010 {
01011         register long32 tmp;
01012         tmp =  (dw & (ulong32)0x000000FF);
01013         tmp = ((dw & (ulong32)0x0000FF00) >> 0x08) | (tmp << (ulong32)0x08);
01014         tmp = ((dw & (ulong32)0x00FF0000) >> 0x10) | (tmp << (ulong32)0x08);
01015         tmp = ((dw & (ulong32)0xFF000000) >> 0x18) | (tmp << (ulong32)0x08);
01016         return(tmp);
01017 }
01018 
01019 #endif
01020 
01021 
01022 /*-------------------------------------------------------------------------*/
01046 /*--------------------------------------------------------------------------*/
01047 
01048 void
01049 convert_fits_to_local(
01050         pixelvalue      *       p_dest,
01051         BYTE            *       p_source,
01052         ulong32                             nbpix,
01053         int                             pixel_type,
01054         double                  b_scale,
01055         double                  b_zero)
01056 {
01057         int                     i ;
01058         short16         shortintpix ;
01059         long32          longintpix ;
01060         float           floatpix ;
01061         double          doublepix ;
01062         BYTE            XLpix[8] ;
01063 
01064         /* Error handling : test entries */
01065 
01066         if (p_dest == NULL || p_source == NULL) return ;
01067 
01068         if ((nbpix==0)||(nbpix > MAX_COLUMN_NUMBER * MAX_LINE_NUMBER)) {
01069                 e_error("wrong number of pixels: %ld", nbpix) ;
01070                 return ;
01071         }
01072 
01073         if (b_scale < 1e-50) {
01074                 e_warning("wrong value for BSCALE: %g", b_scale) ;
01075         }
01076 
01077         switch(pixel_type) {
01078                 case BPP_8_UNSIGNED:
01079                 /* No swapping needed : we are dealing with bytes       */
01080                 for (i=0 ; i<nbpix ; i++)
01081                         *p_dest++ = (pixelvalue)(b_scale * (double)(*p_source++) + b_zero);
01082                 break ;
01083 
01084                 case BPP_16_SIGNED:
01085                 /* Swapping needed for 16 bit data if machine is not MOTOROLA */
01086 #if (BYTE_ORDERING_SCHEME == INTEL)
01087                 /* INTEL like based processors  */
01088                 for (i=0 ; i<nbpix ; i++) {
01089                         memcpy(&shortintpix, p_source, 2) ;
01090                         p_source += 2 ;
01091                         shortintpix = swap_2_bytes(shortintpix) ;
01092                         *p_dest++ = (pixelvalue)(b_scale * (double)shortintpix + b_zero);
01093                 }
01094 #else
01095                 /* MOTOROLA like based processors       */
01096                 for (i=0 ; i<nbpix ; i++) {
01097                         memcpy(&shortintpix, p_source, 2) ;
01098                         p_source += 2 ;
01099                         *p_dest++ = (pixelvalue)(b_scale * (double)shortintpix + b_zero);
01100                 }
01101 #endif
01102                 break ;
01103 
01104                 case BPP_32_SIGNED:
01105                 /* Swapping needed for 32 bit data if machine is not MOTOROLA */
01106 #if (BYTE_ORDERING_SCHEME == INTEL)
01107                 /* INTEL like based processors  */
01108                 for (i=0 ; i<nbpix ; i++) {
01109                         memcpy(&longintpix, p_source, 4) ;
01110                         p_source += 4 ;
01111                         longintpix = swap_4_bytes(longintpix) ;
01112                         *p_dest++ = (pixelvalue)(b_scale * (double)longintpix + b_zero) ;
01113                 }
01114 #else
01115                 /* MOTOROLA like based processors       */
01116                 for (i=0 ; i<nbpix ; i++) {
01117                         memcpy(&longintpix, p_source, 4) ;
01118                         p_source += 4 ;
01119                         *p_dest++ = (pixelvalue)(b_scale * (double)longintpix + b_zero) ;
01120                 }
01121 #endif
01122 
01123                 break ;
01124 
01125                 case BPP_IEEE_FLOAT:
01126                 /* Swapping needed for 32 bit data if machine is not MOTOROLA */
01127 #if (BYTE_ORDERING_SCHEME == INTEL)
01128                 /* INTEL like based processors  */
01129                 for (i=0 ; i<nbpix ; i++) {
01130                         memcpy(&longintpix, p_source, 4) ;
01131                         p_source += 4 ;
01132                         longintpix = swap_4_bytes(longintpix) ;
01133                         memcpy(&floatpix, &longintpix, 4) ;
01134                         *p_dest++ = (pixelvalue)(b_scale * (double)floatpix + b_zero) ;
01135                 }
01136 #else
01137                 /* MOTOROLA like based processors       */
01138                 for (i=0 ; i<nbpix ; i++) {
01139                         memcpy(&longintpix, p_source, 4) ;
01140                         p_source += 4 ;
01141                         memcpy(&floatpix, &longintpix, 4) ;
01142                         *p_dest++ = (pixelvalue)(b_scale * (double)floatpix + b_zero) ;
01143                 }
01144 #endif
01145                 break ;
01146 
01147 
01148                 case BPP_IEEE_DOUBLE:
01149                 /*
01150                  * Double data is converted to internal format for
01151                  * storage: beware of loss of precision !!
01152                  *
01153                  * Still: double data has to be swapped
01154                  */
01155                 if (sizeof(pixelvalue)<sizeof(double)) {
01156                         e_warning("loss of precision in input") ;
01157                 }
01158 #if (BYTE_ORDERING_SCHEME == INTEL)
01159                 /* INTEL like based processors
01160                  * endian garbaging does the following on 64 bits:
01161                  *
01162                  * 1234567890ABCDEF becomes
01163                  * EFCDAB9078563412
01164                  */
01165                 for (i=0 ; i<nbpix ; i++) {
01166                         /* Copy input bytes into 8 bytes, no loop to optimize */
01167                         XLpix[0] = *p_source++ ;
01168                         XLpix[1] = *p_source++ ;
01169                         XLpix[2] = *p_source++ ;
01170                         XLpix[3] = *p_source++ ;
01171                         XLpix[4] = *p_source++ ;
01172                         XLpix[5] = *p_source++ ;
01173                         XLpix[6] = *p_source++ ;
01174                         XLpix[7] = *p_source++ ;
01175                         /* swap_8_bytes(XLpix) ; */
01176                         swap_bytes(XLpix, 8) ;
01177                         doublepix = *((double*)XLpix) ;
01178                         *p_dest++ =
01179                                 (pixelvalue)((double)b_scale*doublepix+(double)b_zero);
01180                 }
01181 #else
01182                 /* MOTOROLA like based processors       */
01183                 for (i=0 ; i<nbpix ; i++) {
01184                         /* Copy input bytes into 8 bytes, no loop to optimize */
01185                         XLpix[0] = *p_source++ ;
01186                         XLpix[1] = *p_source++ ;
01187                         XLpix[2] = *p_source++ ;
01188                         XLpix[3] = *p_source++ ;
01189                         XLpix[4] = *p_source++ ;
01190                         XLpix[5] = *p_source++ ;
01191                         XLpix[6] = *p_source++ ;
01192                         XLpix[7] = *p_source++ ;
01193                         doublepix = *((double*)XLpix) ;
01194                         *p_dest++ =
01195                                 (pixelvalue)((double)b_scale*doublepix+(double)b_zero);
01196                 }
01197 #endif
01198                 break ;
01199 
01200                 default:
01201                 e_error("unrecognized pixel type in input: aborting conversion") ;
01202                 break ;
01203         }
01204 }
01205 
01206 
01207 /*-------------------------------------------------------------------------*/
01228 /*--------------------------------------------------------------------------*/
01229 
01230 
01231 OneCube * list_make_cube(OneImage ** list, int np)
01232 {
01233         OneCube * cu ;
01234         int               i ;
01235 
01236         /* Check entries */
01237         if (list==NULL || np<1) return NULL ;
01238 
01239         /* Check image pointers and image sizes */
01240         for (i=0 ; i<np ; i++) {
01241                 if (list[i]==NULL) {
01242                         e_error("NULL image pointer in list: aborting cube creation");
01243                         return NULL ;
01244                 }
01245                 if ((list[i]->lx != list[0]->lx) ||
01246                         (list[i]->ly != list[0]->ly)) {
01247                         e_error("images have different sizes: aborting cube creation");
01248                         return NULL ;
01249                 }
01250         }
01251 
01252         /* Create new cube, fill it up with image pointers */
01253         cu = new_cube(list[0]->lx, list[0]->ly, np);
01254         memcpy(cu->plane, list, np * sizeof(OneImage*));
01255 
01256         return cu ;
01257 }
01258 /*-------------------------------------------------------------------------*/
01273 /*--------------------------------------------------------------------------*/
01274 
01275 OneCube *
01276 new_cube(
01277     int     lx,
01278     int     ly,
01279     int     n_im
01280 )
01281 {
01282     OneCube *   new_cube ;
01283     int         i ;
01284 
01285     if ((lx<=0)||(lx>MAX_COLUMN_NUMBER)||(ly<=0)||(ly>MAX_LINE_NUMBER) ||
01286         (n_im<=0) || (n_im>MAX_IMAGE_NUMBER)) {
01287         e_error("error in requested cube size: [%d x %d x %d]", lx, ly, n_im) ;
01288         return NULL ;
01289     }
01290 
01291     /* Try to get allocation for a new structure    */
01292     new_cube = malloc(sizeof(OneCube)) ;
01293         new_cube->plane = malloc(n_im * sizeof(OneImage*)) ;
01294     for (i=0 ; i<n_im ; i++) {
01295         new_cube->plane[i] = NULL ;
01296         }
01297     new_cube->type = 0;
01298     new_cube->lx = lx ;
01299     new_cube->ly = ly ;
01300     new_cube->np = n_im ;
01301     new_cube->nbpix = lx * ly * n_im ;
01302     new_cube->history = NULL ;
01303     new_cube->n_comments = 0 ;
01304     new_cube->orig_ptype = BPP_DEFAULT ;
01305         new_cube->filename = NULL ;
01306     new_cube->data=NULL;
01307     new_cube->planes=NULL;
01308 
01309     return new_cube ;
01310 }
01311 
01312 
01313 /*-------------------------------------------------------------------------*/
01325 /*--------------------------------------------------------------------------*/
01326 
01327 #define PIX_SWAP(a,b) { pixelvalue temp=(a);(a)=(b);(b)=temp; }
01328 #define PIX_STACK_SIZE 50
01329 
01330 void pixel_qsort(pixelvalue *pix_arr, int npix)
01331 {
01332     int         i,
01333                 ir,
01334                 j,
01335                 k,
01336                 l;
01337     int        i_stack[PIX_STACK_SIZE*sizeof(pixelvalue)] ;
01338     int         j_stack ;
01339     pixelvalue  a ;
01340 
01341     ir = npix ;
01342     l = 1 ;
01343     j_stack = 0 ;
01344     for (;;) {
01345         if (ir-l < 7) {
01346             for (j=l+1 ; j<=ir ; j++) {
01347                 a = pix_arr[j-1];
01348                 for (i=j-1 ; i>=1 ; i--) {
01349                     if (pix_arr[i-1] <= a) break;
01350                     pix_arr[i] = pix_arr[i-1];
01351                 }
01352                 pix_arr[i] = a;
01353             }
01354             if (j_stack == 0) break;
01355             ir = i_stack[j_stack-- -1];
01356             l  = i_stack[j_stack-- -1];
01357         } else {
01358             k = (l+ir) >> 1;
01359             PIX_SWAP(pix_arr[k-1], pix_arr[l])
01360             if (pix_arr[l] > pix_arr[ir-1]) {
01361                 PIX_SWAP(pix_arr[l], pix_arr[ir-1])
01362             }
01363             if (pix_arr[l-1] > pix_arr[ir-1]) {
01364                 PIX_SWAP(pix_arr[l-1], pix_arr[ir-1])
01365             }
01366             if (pix_arr[l] > pix_arr[l-1]) {
01367                 PIX_SWAP(pix_arr[l], pix_arr[l-1])
01368             }
01369            i = l+1;
01370             j = ir;
01371             a = pix_arr[l-1];
01372             for (;;) {
01373                 do i++; while (pix_arr[i-1] < a);
01374                 do j--; while (pix_arr[j-1] > a);
01375                 if (j < i) break;
01376                 PIX_SWAP(pix_arr[i-1], pix_arr[j-1]);
01377             }
01378             pix_arr[l-1] = pix_arr[j-1];
01379             pix_arr[j-1] = a;
01380             j_stack += 2;
01381             if (j_stack > PIX_STACK_SIZE) {
01382                 e_error("stack too small in pixel_qsort: aborting");
01383                 exit(-2001) ;
01384             }
01385             if (ir-i+1 >= j-l) {
01386                 i_stack[j_stack-1] = ir;
01387                 i_stack[j_stack-2] = i;
01388                 ir = j-1;
01389             } else {
01390                 i_stack[j_stack-1] = j-1;
01391                 i_stack[j_stack-2] = l;
01392                 l = i;
01393             }
01394         }
01395     }
01396 }
01397 #undef PIX_STACK_SIZE
01398 #undef PIX_SWAP
01399 
01400 /*-------------------------------------------------------------------------*/
01410 /*--------------------------------------------------------------------------*/
01411 void
01412 destroy_image(OneImage * d)
01413 {
01414     if (d == NULL) return ;
01415         if (d->reg.active==0) {
01416                 /* Image is in memory: de-allocate properly */
01417         if (d->data != NULL)
01418             free(d->data) ;
01419         } else {
01420                 /*
01421                  * Image comes from an mmap'ed file, register out the image
01422                  * from the reg_alloc table. If the count has reached zero,
01423                  * the file can now be unmapped.
01424                  */
01425                 if (reg_alloc_remove(&(d->reg))==0) {
01426                         mmap_close(d->reg.mm);
01427                 }
01428         }
01429         free(d) ;
01430 }
01431 /*-------------------------------------------------------------------------*/
01445 /*--------------------------------------------------------------------------*/
01446 
01447 void
01448 destroy_cube(
01449     OneCube * d
01450 )
01451 {
01452     int     i ;
01453 
01454     if (d == NULL) return ;
01455 
01456         /* Destroy all attached images */
01457         if (d->type == 1){
01458         if (d->data != NULL)
01459                 free(d->data) ;
01460         if (d->planes != NULL)
01461                 free(d->planes) ;
01462         }
01463         else if (d->np>0) {
01464                 for (i=0 ; i<d->np ; i++) {
01465                         if (d->plane[i]!=NULL) {
01466                                 destroy_image(d->plane[i]);
01467                         }
01468                 }
01469         }
01470 
01471         /* free general plane pointer */
01472         if (d->plane != NULL)
01473                 free(d->plane) ;
01474 
01475     /* free history zone    */
01476         if (d->n_comments>0) {
01477                 for (i=0 ; i<d->n_comments ; i++) {
01478                         free(d->history[i]);
01479                 }
01480                 free(d->history);
01481         }
01482     /* free structure itself    */
01483     free(d) ;
01484 
01485 }
01486 
01487 /*-------------------------------------------------------------------------*/
01507 /*--------------------------------------------------------------------------*/
01508 
01509 OneImage *
01510 average_with_rejection(
01511         OneCube *       incube,
01512         double          lo_reject,
01513         double          hi_reject
01514 )
01515 {
01516         OneImage        *       avg ;
01517         int                             pos ;
01518         int                             i, j ;
01519         pixelvalue              rejavg ;
01520 
01521         /* Error handling: test entries */
01522         if (incube==NULL) return NULL ;
01523         if (incube->np<3) {
01524                 /*e_warning("no average rejection with < 3 planes - simple averaging") ;*/
01525                 avg = average_cube_to_image(incube) ;
01526                 if (avg == NULL) {
01527                         e_error("cannot compute the simple average") ;
01528                         return NULL ;
01529                 }
01530         } else {
01531 
01532                 if ((lo_reject + hi_reject) > 0.9) {
01533                         e_error("illegal rejection thresholds: [%f] and [%f]",
01534                                   lo_reject, hi_reject) ;
01535                         e_error("threshold sum should not be over 90%: aborting average") ;
01536                         return NULL ;
01537                 }
01538 
01539                 avg = new_image(incube->lx, incube->ly) ;
01540                 for (j=0 ; j<incube->ly ; j++) {
01541                     compute_status("averaging with rejection", j, incube->ly, 1) ;
01542                         for (i=0 ; i<incube->lx ; i++) {
01543                                 pos = i + j * incube->lx ;
01544                                 rejavg = average_reject_on_time_line(incube,
01545 
01546 pos,
01547 
01548 lo_reject,
01549 
01550 hi_reject) ;
01551                                 avg->data[pos] = rejavg ;
01552                         }
01553                 }
01554         }
01555         return(avg) ;
01556 }
01557 
01558 
01559 
01560 /*-------------------------------------------------------------------------*/
01575 /*--------------------------------------------------------------------------*/
01576 
01577 
01578 void
01579 e_warning(const char *fmt, ...)
01580 {
01581         char   *msg ;
01582         va_list ap ;
01583         char    logmsg[1024] ;
01584 
01585         fprintf(stderr, "*** ") ;
01586         msg = strdup(fmt) ;
01587         chop_r(msg) ;
01588 
01589         va_start(ap, fmt) ;
01590         vfprintf(stderr, msg, ap) ;
01591         vsprintf(logmsg, msg, ap) ;
01592         va_end(ap) ;
01593 
01594         fprintf(stderr, " ***\n") ;
01595         fflush(stderr) ;
01596 
01597         free(msg) ;
01598 
01599         if (eclipse_logfile) {
01600                 e_logfile("warning", logmsg);
01601         }
01602         return ;
01603 }
01604 
01605 /*-------------------------------------------------------------------------*/
01621 /*--------------------------------------------------------------------------*/
01622 
01623 
01624 void
01625 e_error(const char *fmt, ...)
01626 {
01627         char  * msg ;
01628         va_list ap ;
01629         char    logmsg[1024] ;
01630 
01631         fprintf(stderr, "error: ") ;
01632         msg = strdup(fmt) ;
01633         chop_r(msg) ;
01634 
01635         va_start(ap, fmt) ;
01636         vfprintf(stderr, msg, ap) ;
01637         vsprintf(logmsg, msg, ap) ;
01638         va_end(ap) ;
01639 
01640         free(msg) ;
01641         fprintf(stderr, "\n") ;
01642         fflush(stderr) ;
01643 
01644         if (eclipse_logfile) {
01645                 e_logfile("error", logmsg);
01646         }
01647 }
01648 
01649 /*-------------------------------------------------------------------------*/
01667 /*--------------------------------------------------------------------------*/
01668 
01669 void
01670 e_comment(int level, const char * fmt, ...)
01671 {
01672         int     i ;
01673         char  * msg ;
01674         va_list ap ;
01675         char    logmsg[1024] ;
01676 
01677         if (eclipse_verbose) {
01678                 if (level != 0)
01679                         for (i=0 ; i<level ; i++)
01680                                 fprintf(stderr, "\t") ;
01681 
01682                 msg = strdup(fmt) ;
01683                 chop_r(msg) ;
01684                 va_start(ap, fmt);
01685                 vfprintf(stderr, msg, ap) ;
01686                 vsprintf(logmsg, msg, ap) ;
01687                 va_end(ap) ;
01688 
01689                 free(msg) ;
01690                 fprintf(stderr,"\n") ;
01691                 fflush(stderr) ;
01692                 if (eclipse_logfile) {
01693                         e_logfile("info", logmsg);
01694                 }
01695         }
01696 }
01697 
01698 /*-------------------------------------------------------------------------*/
01744 /*--------------------------------------------------------------------------*/
01745 
01746 void
01747 compute_status(char *Msg, int done, int total, int level)
01748 {
01749         int     i ;
01750         char    logmsg[1024];
01751 
01752         if (eclipse_verbose) {
01753                 done++ ;
01754                 printf("\r");
01755                 if (level) {
01756                        for (i=0 ; i<level ; i++) {
01757                                 printf("\t");
01758                         }
01759                 }
01760                 printf("%s: %d out of %d ", Msg, done, total);
01761                 fflush(stdout) ;
01762 
01763                 if (done >= total)
01764                         fprintf(stdout, "\n") ;
01765         }
01766         if (done==total && eclipse_verbose) {
01767                 if (eclipse_logfile) {
01768                         sprintf(logmsg, "%s: completed", Msg);
01769                         e_logfile("info", logmsg);
01770                 }
01771         }
01772         return ;
01773 }
01774 
01775 /*-------------------------------------------------------------------------*/
01797 /*--------------------------------------------------------------------------*/
01798 
01799 char * fits_pretty_string(char * s)
01800 {
01801     static char     pretty[81] ;
01802     int             i,j ;
01803 
01804         /* bulletproof */
01805         if (s==NULL) return NULL ;
01806 
01807     pretty[0] = (char)0 ;
01808     if (s[0]!='\'') return s ;
01809 
01810     /* skip first quote */
01811     i=1 ;
01812     j=0 ;
01813     /* trim left-side blanks */
01814     while (s[i]==' ') {
01815         if (i==(int)strlen(s)) break ;
01816         i++ ;
01817     }
01818     if (i>=(int)(strlen(s)-1)) return pretty ;
01819     /* copy string, changing double quotes to single ones */
01820     while (i<(int)strlen(s)) {
01821         if (s[i]=='\'') {
01822             i++ ;
01823         }
01824         pretty[j]=s[i];
01825         i++ ;
01826         j++ ;
01827     }
01828     /* NULL-terminate the pretty string */
01829     pretty[j+1]=(char)0;
01830     /* trim right-side blanks */
01831     j = (int)strlen(pretty)-1;
01832     while (pretty[j]==' ') j-- ;
01833     pretty[j+1]=(char)0;
01834     return pretty;
01835 }
01836 
01837 /*-------------------------------------------------------------------------*/
01848 /*--------------------------------------------------------------------------*/
01849 
01850 OneImage        *
01851 average_cube_to_image(
01852         OneCube *       incube
01853 )
01854 {
01855         OneImage        *       sum_image ;
01856         int                             i ;
01857         pixelvalue              inv ;
01858 
01859         if (incube==NULL) return NULL ;
01860         e_comment(1, "averaging cube to one image") ;
01861         sum_image = new_image(incube->lx, incube->ly) ;
01862         if (sum_image==NULL) return NULL ;
01863         /* Loop on all planes */
01864         for (i=0 ; i<incube->np ; i++) {
01865                 compute_status("linear averaging", i, incube->np, 2) ;
01866                 add_img_local(sum_image, incube->plane[i]) ;
01867         }
01868         /* average on number of planes */
01869         inv = (pixelvalue)(1.0 / (double)incube->np) ;
01870         for (i=0 ; i<sum_image->nbpix ; i++) {
01871                 sum_image->data[i] *= inv ;
01872         }
01873         return(sum_image) ;
01874 }
01875 
01876 /*-------------------------------------------------------------------------*/
01897 /*--------------------------------------------------------------------------*/
01898 
01899 pixelvalue
01900 average_reject_on_time_line(
01901         OneCube         *       in_cube,
01902         ulong32                             pos,
01903         double                  lo_rej,
01904         double                  hi_rej)
01905 {
01906         int                             plane ;
01907         pixelvalue              timeline[2048] ;
01908         double                  acc_val ;
01909         int                             lo_n,
01910                                         hi_n ;
01911         int                             i ;
01912         int                             nv ;
01913 
01914         if (in_cube == NULL) return (pixelvalue)-1 ;
01915         if (in_cube->np < 3) {
01916                 e_error("median average has no meaning with less than 3 planes");
01917                 return (pixelvalue)-1 ;
01918         }
01919 
01920         lo_n = (int)(in_cube->np * lo_rej +0.5);
01921         hi_n = (int)(in_cube->np * hi_rej +0.5);
01922         if (lo_n + hi_n >= in_cube->np) {
01923                 return 0.0 ;
01924         }
01925 
01926 
01927         /*
01928          * Get all the pixel values along timeline
01929          */
01930 
01931         for (plane=0 ; plane<in_cube->np ; plane++) {
01932                 timeline[plane] = in_cube->plane[plane]->data[pos] ;
01933         }
01934        /* Now sort out the timeline for this pixel     */
01935         pixel_qsort(timeline, in_cube->np) ;
01936         acc_val = 0.0 ;
01937 
01938         /*
01939          * Get only the middle values, rejecting the lower and upper pixel
01940          * proportion as requested
01941          */
01942 
01943         nv = 0 ;
01944         for (i=lo_n ; i<(in_cube->np - hi_n) ; i++) {
01945                 acc_val += (double)timeline[i] ;
01946                 nv++ ;
01947         }
01948         acc_val /= (double)nv ;
01949         return (pixelvalue)acc_val ;
01950 }
01951 
01952 /*-------------------------------------------------------------------------*/
01969 /*--------------------------------------------------------------------------*/
01970 
01971 void e_logfile(char * type, char * msg)
01972 {
01973         static int      logfile_started = 0 ;
01974         FILE    *       log ;
01975 
01976         if ((log = fopen(eclipse_logfilename, "a")) == NULL) return ;
01977 
01978         if (logfile_started==0) {
01979                 e_logfile_start();
01980                 logfile_started = 1 ;
01981         }
01982 
01983         fprintf(log, "[%s]: %s\n", type, msg);
01984         fclose(log);
01985         return ;
01986 }
01987 
01988 
01989 /*-------------------------------------------------------------------------*/
02005 /*--------------------------------------------------------------------------*/
02006 
02007 void e_logfile_start(void)
02008 {
02009         FILE    *       log ;
02010         char            hostname[HOSTNAMESZ] ;
02011     char        time_str[TIMESTRSZ] ;
02012     time_t      local_t ;
02013 
02014         /*
02015          * Open LOGFILE in append
02016          */
02017 
02018         if ((log = fopen(eclipse_logfilename, "a"))==NULL) {
02019                 fprintf(stderr, "error: cannot open log file\n");
02020                 return ;
02021         }
02022 
02023         /*
02024          * Get hostname, username, PID, date/time
02025          */
02026     if (gethostname(hostname, (size_t)HOSTNAMESZ) == -1) {
02027         sprintf(hostname, "localhost") ;
02028     }
02029     if (time(&local_t) == (time_t)-1) {
02030         sprintf(time_str, "unknown date and time") ;
02031     } else {
02032         sprintf(time_str, "%s", ctime(&local_t)) ;
02033     }
02034 
02035         fprintf(log, "\n\n\n");
02036         fprintf(log, "->---------------------------------------------------\n");
02037         fprintf(log, "-> start logging\n");
02038         fprintf(log, "-> [%s@%s] PID [%ld]\n",
02039                         get_login_name(),
02040                         hostname,
02041                         (long)getpid());
02042         fprintf(log, "-> %s\n", time_str);
02043         fprintf(log, "\n\n") ;
02044 
02045         fclose(log);
02046        return ;
02047 }
02048 
02049 /*-------------------------------------------------------------------------*/
02064 /*--------------------------------------------------------------------------*/
02065 
02066 void e_logfile_stop(void)
02067 {
02068         FILE    *       log;
02069     char        time_str[TIMESTRSZ] ;
02070     time_t      local_t ;
02071 
02072     /*
02073      * Open LOGFILE in append
02074      */
02075 
02076     if ((log = fopen(eclipse_logfilename, "a"))==NULL) {
02077         fprintf(stderr, "error: cannot open log file\n");
02078         return ;
02079     }
02080 
02081     if (time(&local_t) == (time_t)-1) {
02082         sprintf(time_str, "unknown date and time") ;
02083     } else {
02084         sprintf(time_str, "%s", ctime(&local_t)) ;
02085     }
02086 
02087         fprintf(log, "\n");
02088         fprintf(log, "-> stop logging\n");
02089         fprintf(log, "-> %s\n", time_str);
02090         fclose(log);
02091         return ;
02092 }
02093 
02094 
02095 /*-------------------------------------------------------------------------*/
02107 /*--------------------------------------------------------------------------*/
02108 
02109 
02110 int
02111 add_img_local(
02112         OneImage * im1,
02113         OneImage * im2
02114 )
02115 {
02116         register pixelvalue * p1,
02117                                                 * p2 ;
02118         register int              i ;
02119 
02120         p1 = im1->data ;
02121         p2 = im2->data ;
02122         for (i=0 ; i<im1->nbpix ; i++) {
02123                 *p1++ += *p2++ ;
02124         }
02125         return 0 ;
02126 }
02127 
02128 
02129 
02130 
02131 /*-------------------------------------------------------------------------*/
02148 /*--------------------------------------------------------------------------*/
02149 
02150 void memory_status(void)
02151 {
02152         malloc_cell     *curr_cell ;
02153 
02154         if ((total_allocated_RAM!=0) ||
02155                 (total_allocated_SWAP!=0) ||
02156                 (nswapfiles>0) ||
02157                 (nalloc_ptrs>0)) {
02158                 (void)fprintf(stderr, "**** MEMORY STATUS\n") ;
02159                 (void)fprintf(stderr, "memory in use  : %d bytes\n",
02160                                 total_allocated_RAM + total_allocated_SWAP) ;
02161                 (void)fprintf(stderr, "RAM            : %d bytes\n",
02162                                 total_allocated_RAM);
02163                 (void)fprintf(stderr, "SWAP           : %d bytes\n",
02164                                 total_allocated_SWAP);
02165                 (void)fprintf(stderr, "open swapfiles : %d\n", nswapfiles) ;
02166 
02167                 curr_cell = malloc_table ;
02168                 while (curr_cell != NULL) {
02169 
02170                         if (curr_cell->flavour == ALLOC_CELL_RAM) {
02171                                 (void)fprintf(stderr, "[RAM] ") ;
02172                         } else {
02173                                 (void)fprintf(stderr, "[SWP] ") ;
02174                         }
02175 
02176                         (void)fprintf(stderr, "po=%p size=%ld allocated in %s (%d) ",
02177                                         curr_cell->pointer,
02178                                         (long)curr_cell->size,
02179                                         curr_cell->filename,
02180                                         curr_cell->lineno) ;
02181 
02182                         if (curr_cell->flavour == ALLOC_CELL_SWAP) {
02183                                 (void)fprintf(stderr, "in swapfile %s",
02184                                                 curr_cell->swapfilename) ;
02185                         }
02186                         (void)fprintf(stderr, "\n") ;
02187                         (void)fflush(stderr) ;
02188 
02189                         curr_cell = curr_cell->next ;
02190                 }
02191         }
02192         return ;
02193 }
02194 
02195 /*-------------------------------------------------------------------------*/
02207 /*--------------------------------------------------------------------------*/
02208 
02209 OneImage *
02210 sub_image(
02211         OneImage *image1,
02212         OneImage *image2)
02213 {
02214     OneImage    *image_out ;
02215     size_t              i;
02216         register
02217         pixelvalue              *p1, *p2, *outp ;
02218 
02219         if (image1==NULL || image2==NULL) return NULL ;
02220 
02221     /* Input data images shall have the same sizes */
02222     if (    (image1->lx != image2->lx) ||
02223             (image1->ly != image2->ly)    ) {
02224                 e_error("image1 is [%d x %d] image2 is [%d x %d]",
02225                                 image1->lx, image1->ly, image2->lx, image2->ly) ;
02226                 e_error("cannot subtract images of different sizes\n") ;
02227                 return NULL ;
02228     }
02229 
02230     /* Allocations  */
02231     image_out = new_image(image1->lx, image1->ly) ;
02232 
02233     /* Subtract two images */
02234         p1 = image1->data ;
02235         p2 = image2->data ;
02236         outp = image_out->data ;
02237 
02238     for (i=0 ; i<image_out->nbpix ; i++)
02239                 *outp++ = *p1++ - *p2++ ;
02240 
02241     return image_out ;
02242 }
02243 
02244 /*-------------------------------------------------------------------------*/
02255 /*--------------------------------------------------------------------------*/
02256 
02257 OneImage *
02258 copy_image(OneImage * src_img)
02259 {
02260     OneImage * dest_img ;
02261 
02262     if (src_img==NULL) return NULL ;
02263     dest_img = new_image(src_img->lx, src_img->ly) ;
02264     memcpy(dest_img->data,
02265            src_img->data,
02266            (size_t)src_img->nbpix * sizeof(pixelvalue)) ;
02267     return dest_img ;
02268 }
02269 
02270 /*-------------------------------------------------------------------------*/
02282 /*--------------------------------------------------------------------------*/
02283 
02284 char *
02285 get_tempfile_name(void)
02286 {
02287         static char tmpfilename[TMPFILENAME_SZ] ;
02288         file_reg++ ;
02289         sprintf(tmpfilename, "%s/vmswap_%05ld_%05x",
02290                           emem_tmpdirname,
02291                           (long)getpid(),
02292                           file_reg) ;
02293         return tmpfilename ;
02294 }
02295 
02296 
02297 
02298 /*-------------------------------------------------------------------------*/
02320 /*--------------------------------------------------------------------------*/
02321 
02322 static void *
02323 swap_malloc_x(
02324         size_t          size,
02325         char    *       filename,
02326         int                     lineno
02327 )
02328 {
02329         void            *       pointer ;
02330         malloc_cell     *       new_cell ;
02331         malloc_cell     *       curr_cell ;
02332 
02333         char            *       fname ;
02334         int                             swapfiled ;
02335         char            *       wbuf ;
02336         int                             nbufs ;
02337         int                             i ;
02338 
02339         if (size+total_allocated_SWAP > emem_vm_limit*ONE_MEG) {
02340                 e_error("fatal: swap space overflow --> killing process\n") ;
02341                 exit(-2000) ;
02342         }
02343 
02344         wbuf = calloc(emem_vm_pagesize,1) ;
02345         if (wbuf == NULL) {
02346                 e_error("fatal: cannot allocate internal buffer\n") ;
02347                 e_error("exiting program now\n") ;
02348                 exit(-2000) ;
02349         }
02350 
02351         /* create swap file with rights: rwxrwxrwx */
02352         fname = get_tempfile_name() ;
02353         swapfiled = open(fname, O_RDWR | O_CREAT) ;
02354         fchmod(swapfiled, S_IRWXU | S_IRWXG | S_IRWXO) ;
02355 
02356     if (size % emem_vm_pagesize == 0) {
02357         nbufs = size / emem_vm_pagesize ;
02358     } else {
02359         nbufs = 1 + (size / emem_vm_pagesize) ;
02360     }
02361     for (i=0 ; i<nbufs ; i++) {
02362         if (write(swapfiled, wbuf, emem_vm_pagesize) == -1) {
02363                         perror("write") ;
02364                         e_error("eclipse\n");
02365                         e_error("fatal: cannot write to swapfile: [%s] full?\n", emem_tmpdirname);
02366                         e_error("exiting program now\n") ;
02367                         free(wbuf) ;
02368                         close(swapfiled) ;
02369                         remove(fname) ;
02370                         free(fname) ;
02371                         exit(-2000) ;
02372                 }
02373     }
02374         free(wbuf) ;
02375 
02376     /* mmap() the swap file */
02377     pointer = (void*)mmap(0,
02378                           size,
02379                           PROT_READ | PROT_WRITE,
02380                           MAP_PRIVATE,
02381                           swapfiled,
02382                           0) ;
02383     if (pointer == (void*)-1) {
02384         e_error("mmap failed with errno = %d\n", errno) ;
02385                 e_error("failed to allocate SWAP %ld bytes in %s (%d)\n",
02386                                 size, filename, lineno) ;
02387                 e_error("total SWAP allocated so far : %ld bytes\n",total_allocated_SWAP);
02388                 e_error("%d swap files opened so far in %s\n", nswapfiles, emem_tmpdirname) ;
02389         close(swapfiled) ;
02390                 e_error("exiting program now\n") ;
02391                 exit(-2000) ;
02392     }
02393 #ifdef DEBUG_ALLOC
02394         (void)fprintf(stderr,
02395                                   "-> %p = SWAP malloc(%ld) by %s (%d) in %s/%s\n",
02396                                   pointer, size, filename, lineno, emem_tmpdirname, fname) ;
02397         (void)fflush(stderr) ;
02398 #endif
02399 
02400         total_allocated_SWAP += size ;
02401         nalloc_ptrs ++ ;
02402         nswapfiles++ ;
02403 
02404         new_cell = malloc(sizeof(malloc_cell)) ;
02405         if (new_cell == NULL) {
02406                 e_error("fatal error in memory: life functions terminated\n") ;
02407                 exit(2001) ;
02408         }
02409         new_cell->pointer = pointer ;
02410         new_cell->size = size ;
02411         strncpy(new_cell->filename, filename, SRCFILENAME_SZ-1) ;
02412        new_cell->lineno = lineno ;
02413         new_cell->flavour = ALLOC_CELL_SWAP ;
02414         strncpy(new_cell->swapfilename, fname, TMPFILENAME_SZ-1) ;
02415         new_cell->swapfd = swapfiled ;
02416         new_cell->next = NULL ;
02417         new_cell->prev = NULL ;
02418 
02419         if (malloc_table == NULL) {
02420                 malloc_table = new_cell ;
02421         } else {
02422                 curr_cell = malloc_table ;
02423                 while (curr_cell->next != NULL) {
02424                         curr_cell = curr_cell->next ;
02425                 }
02426                 curr_cell->next = new_cell ;
02427                 new_cell->prev  = curr_cell ;
02428         }
02429         return(pointer) ;
02430 }
02431 
02432 /*-------------------------------------------------------------------------*/
02454 /*--------------------------------------------------------------------------*/
02455 
02456 void *
02457 ext_malloc_x(
02458         size_t          size,
02459         char    *       filename,
02460         int                     lineno
02461 )
02462 {
02463         void    *       pointer ;
02464 
02465         /*
02466          * See if we need to allocate swap for the requested size.
02467          */
02468         if (size > emem_vm_pagesize) {
02469                 if ((size+total_allocated_RAM)>emem_ramlimit*ONE_MEG) {
02470                         pointer = swap_malloc_x(size, filename, lineno) ;
02471                         return pointer ;
02472                 }
02473         }
02474 
02475         /*
02476          * Otherwise allocate in RAM
02477          */
02478         pointer = ram_malloc(size, filename, lineno) ;
02479         return pointer ;
02480 }
02481 
02482 /*-------------------------------------------------------------------------*/
02497 /*--------------------------------------------------------------------------*/
02498 
02499 char * get_login_name(void)
02500 {
02501     static char name[32];
02502     FILE * whoami ;
02503 
02504         name[0]=(char)0 ;
02505     if ((whoami = popen("whoami", "r"))==NULL) return name ;
02506     if ((fgets(name, 31, whoami))==NULL) {
02507                 pclose(whoami);
02508                 return name ;
02509         }
02510     if ((pclose(whoami))==-1) return name ;
02511         if (name[0]!=(char)0) name[strlen(name)-1]=(char)0 ;
02512     return name;
02513 
02514 }
02515 /*-------------------------------------------------------------------------*/
02535 /*--------------------------------------------------------------------------*/
02536 
02537 
02538 static void *
02539 ram_malloc(
02540         size_t          size,
02541         char    *       filename,
02542         int                     lineno
02543 )
02544 {
02545         void            *pointer ;
02546         malloc_cell     *new_cell ;
02547         malloc_cell     *curr_cell ;
02548 
02549         pointer = malloc(size) ;
02550 #ifdef DEBUG_ALLOC
02551         (void)fprintf(stderr,
02552                                   "-> %p = RAM malloc(%ld) in %s (%d)\n",
02553                                   pointer, size, filename, lineno) ;
02554         (void)fflush(stderr) ;
02555 #endif
02556         if (pointer == NULL) {
02557                 e_error("failed to allocate RAM %ld bytes in %s (%d)\n",
02558                                 size, filename, lineno) ;
02559                 e_error("total RAM allocated so far : %ld bytes\n", total_allocated_RAM);
02560                 e_error("exiting program now\n") ;
02561                 exit(-2000) ;
02562         }
02563 
02564         total_allocated_RAM += size ;
02565         nalloc_ptrs ++ ;
02566         new_cell = malloc(sizeof(malloc_cell)) ;
02567         if (new_cell == NULL) {
02568                 e_error("fatal error in memory : life functions terminated\n") ;
02569                 exit(2001) ;
02570         }
02571         new_cell->pointer = pointer ;
02572         new_cell->size = size ;
02573         strncpy(new_cell->filename, filename, SRCFILENAME_SZ-1) ;
02574         new_cell->lineno = lineno ;
02575         new_cell->flavour = ALLOC_CELL_RAM ;
02576         new_cell->next = NULL ;
02577         new_cell->prev = NULL ;
02578 
02579         /*
02580          * Place new memory cell into the structure
02581          */
02582 
02583         if (malloc_table == NULL) {
02584                 malloc_table = new_cell ;
02585         } else {
02586                 curr_cell = malloc_table ;
02587                 while (curr_cell->next != NULL) {
02588                         curr_cell = curr_cell->next ;
02589                 }
02590                 curr_cell->next = new_cell ;
02591                 new_cell->prev  = curr_cell ;
02592         }
02593         return(pointer) ;
02594 }
02595 
02596 
02597 /*-------------------------------------------------------------------------*/
02618 /*--------------------------------------------------------------------------*/
02619 
02620 
02621 static void *
02622 ram_calloc(
02623         size_t          n_elem,
02624         size_t          size,
02625         char    *       filename,
02626         int                     lineno
02627 )
02628 {
02629     void        *pointer ;
02630     malloc_cell  *new_cell ;
02631     malloc_cell  *curr_cell ;
02632 
02633     pointer = (void*)calloc(n_elem, size) ;
02634 #ifdef DEBUG_ALLOC
02635         (void)fprintf(stderr,
02636                                   "-> %p = RAM calloc(%ld,%ld) in %s (%d)\n",
02637                                   pointer, n_elem, size, filename, lineno) ;
02638         (void)fflush(stderr) ;
02639 #endif
02640     if (pointer == NULL) {
02641         e_error("failed to allocate RAM %ld bytes in %s (%d)\n",
02642                                 n_elem * size, filename, lineno) ;
02643                 e_error("total RAM allocated so far: %ld bytes\n", total_allocated_RAM);
02644                 e_error("exiting program now\n") ;
02645                 exit(-2000) ;
02646     }
02647    total_allocated_RAM += n_elem * size ;
02648         nalloc_ptrs ++ ;
02649         new_cell = malloc(sizeof(malloc_cell)) ;
02650         if (new_cell == NULL) {
02651                 e_error("fatal error in memory: life functions terminated\n") ;
02652                 exit(2001) ;
02653         }
02654 
02655         new_cell->pointer = pointer ;
02656         new_cell->size = n_elem * size ;
02657         strncpy(new_cell->filename, filename, SRCFILENAME_SZ-1) ;
02658         new_cell->lineno = lineno ;
02659         new_cell->flavour = ALLOC_CELL_RAM ;
02660         new_cell->next = NULL ;
02661         new_cell->prev = NULL ;
02662 
02663         /*
02664          * Place new memory cell into the structure
02665          */
02666 
02667         if (malloc_table == NULL) {
02668                 malloc_table = new_cell ;
02669         } else {
02670                 curr_cell = malloc_table ;
02671                 while (curr_cell->next != NULL) {
02672                         curr_cell = curr_cell->next ;
02673                 }
02674                 curr_cell->next = new_cell ;
02675                 new_cell->prev  = curr_cell ;
02676         }
02677     return(pointer) ;
02678 }
02679 
02680 
02681 /*-------------------------------------------------------------------------*/
02704 /*--------------------------------------------------------------------------*/
02705 
02706 void *
02707 ext_calloc_x(
02708         size_t          n_elem,
02709         size_t          size,
02710         char    *       filename,
02711         int                     lineno
02712 )
02713 {
02714         size_t          reqsize ;
02715         void    *       pointer ;
02716 
02717         /*
02718          * See if we need to allocate swap for the requested size.
02719          */
02720         reqsize = size * n_elem ;
02721         if (reqsize  > emem_vm_pagesize) {
02722                 if ((reqsize+total_allocated_RAM)>emem_ramlimit*ONE_MEG) {
02723                         pointer = swap_malloc_x(reqsize, filename, lineno) ;
02724                         return pointer ;
02725                 }
02726         }
02727         /*
02728          * Otherwise allocate in RAM
02729          */
02730         pointer = ram_calloc(n_elem, size, filename, lineno) ;
02731         return pointer ;
02732 }
02733 
02734 
02735 void
02736 fill_new_image(
02737     int     size_x,
02738     int     size_y,
02739     OneImage    *       new_image,
02740     pixelvalue  *       data
02741 )
02742 {
02743 
02744     if ((size_x<1) ||
02745         (size_x>MAX_COLUMN_NUMBER) ||
02746         (size_y<1) ||
02747         (size_y>MAX_LINE_NUMBER)) {
02748         e_error("requested size: [%dx%d]", size_x, size_y) ;
02749         e_error("cannot create such size: aborting image creation") ;
02750     }
02751 
02752     new_image->lx = size_x ;
02753     new_image->ly = size_y ;
02754     new_image->nbpix = size_x * size_y ;
02755     new_image->data = data;
02756     new_image->reg.active = 0 ;
02757 }
02758 
02759 /*-------------------------------------------------------------------------*/
02781 /*--------------------------------------------------------------------------*/
02782 
02783 double *
02784 fit_1d_poly(
02785         int                     poly_deg,
02786         dpoint  *       list,
02787         int                     np,
02788         double  *       mse
02789 )
02790 {
02791         int                     i, k ;
02792         Matrix          mA, mB, mX ;
02793         double  *       c ;
02794         double          err ;
02795         double          xp, y ;
02796 
02797         if (np<poly_deg+1) {
02798                 e_error("not enough points") ;
02799                 e_error("cannot fit %dth degree polynomial with %d points",
02800                                 poly_deg, np);
02801                 return NULL;
02802         }
02803 
02804         mA = create_mx(poly_deg+1, np) ;
02805         mB = create_mx(1, np) ;
02806 
02807         for (i=0 ; i<np ; i++) {
02808                 mA->m[i] = 1.0 ;
02809                 for (k=1 ; k<=poly_deg ; k++) {
02810                         mA->m[i+k*np] = ipow(list[i].x, k) ;
02811                 }
02812                 mB->m[i] = list[i].y ;
02813         }
02814 
02815         /*
02816          * Solve XA=B by a least-square solution (aka pseudo-inverse).
02817          */
02818         mX = least_sq_mx(mA,mB) ;
02819         /*
02820          * Delete input matrices
02821          */
02822         close_mx(mA) ;
02823         close_mx(mB) ;
02824         /*
02825          * Examine result
02826          */
02827         if (mX==NULL) {
02828                 e_error("cannot fit: non-invertible matrix") ;
02829                 return NULL ;
02830         }
02831         if (debug_active()>1) {
02832                 fprintf(stdout, "# Least-squares matrix:\n") ;
02833                 print_mx(mX, "mX") ;
02834         }
02835 
02836         c = malloc((poly_deg+1)*sizeof(double)) ;
02837         for (i=0 ; i<(poly_deg+1) ; i++) {
02838                 c[i] = mX->m[i] ;
02839         }
02840         close_mx(mX) ;
02841 
02842         /*
02843          * If requested, compute mean squared error
02844          */
02845         if (mse != NULL) {
02846                 err = 0.00 ;
02847                 for (i=0 ; i<np ; i++) {
02848                         y = c[0] ;
02849                         /*
02850                          * Compute the value obtained through the fit
02851                          */
02852                         for (k=1 ; k<=poly_deg ; k++) {
02853                                 xp = ipow(list[i].x, k) ;
02854                                 y += c[k] * xp ;
02855                         }
02856                         /*
02857                          * Subtract from the true value, square, accumulate
02858                          */
02859                         xp   = ipow(list[i].y - y, 2) ;
02860                         err += xp ;
02861                 }
02862                 /* Average the error term */
02863                 err /= (double)np ;
02864                 *mse = err ;
02865         }
02866         return c ;
02867 }
02868 
02869 /*-------------------------------------------------------------------------*/
02877 /*--------------------------------------------------------------------------*/
02878 
02879 int debug_active(void)
02880 { return eclipse_debug ; }
02881 
02882 /*-------------------------------------------------------------------------*/
02893 /*--------------------------------------------------------------------------*/
02894 
02895 Matrix
02896 create_mx(int nr, int nc)
02897 {
02898         Matrix b;
02899         b = (Matrix)calloc(1,sizeof(matrix));
02900         b->m = (double*)calloc(nr*nc,sizeof(double));
02901         b->nr= nr;
02902         b->nc= nc;
02903         return b;
02904 }
02905 
02906 
02907 /*-------------------------------------------------------------------------*/
02917 /*--------------------------------------------------------------------------*/
02918 
02919 Matrix
02920 copy_mx(Matrix a)
02921 {
02922         Matrix b = create_mx(a->nr,a->nc);
02923         if (b!=NULL) {
02924                 register int s = a->nr*a->nc;
02925                 register double *mm = b->m+s;
02926                 register double *am = a->m+s;
02927                 while (s--) *--mm = *--am;
02928         }
02929         return b;
02930 }
02931 
02932 
02933 /*-------------------------------------------------------------------------*/
02943 /*--------------------------------------------------------------------------*/
02944 
02945 void
02946 close_mx(Matrix a)
02947 {
02948         if (a==NULL) return ;
02949         if (a->m != NULL)
02950                 free(a->m);
02951         free(a);
02952         return;
02953 }
02954 
02955 
02956 
02957 /*-------------------------------------------------------------------------*/
02968 /*--------------------------------------------------------------------------*/
02969 
02970 Matrix
02971 mul_mx(Matrix a, Matrix b)
02972 {
02973         Matrix c, d;
02974         int n1=a->nr, n2=a->nc, n3=b->nc;
02975         register double *a0;
02976         register double *c0;
02977         register double *d0;
02978         register int i,j,k;
02979 
02980         if(n2!=b->nr) return NULL;
02981         c = create_mx(n1,n3);
02982         d = transp_mx(b);
02983 
02984         for (i=0,c0=c->m;i<n1;i++)
02985                 for (j=0,d0=d->m;j<n3;j++,c0++)
02986                         for (k=0,*c0=0,a0=a->m+i*n2;k<n2;k++)
02987                                 *c0 += *a0++ * *d0++;
02988         close_mx(d);
02989         return c;
02990 }
02991 
02992 /*-------------------------------------------------------------------------*/
03004 /*--------------------------------------------------------------------------*/
03005 
03006 Matrix
03007 invert_mx(Matrix aa)
03008 {
03009         Matrix bb;
03010         int test=1;
03011 
03012         if(aa->nr!=aa->nc) return NULL;
03013         bb = create_mx(aa->nr,aa->nc);
03014 
03015         if(aa->nr==1)
03016         {
03017                 double det;
03018                 register double ted;
03019                 det= *(aa->m);
03020                 if(dtiny(det)) test=0;
03021                 ted=1./det;
03022                 *(bb->m)=ted;
03023         }
03024         else if(aa->nr==2)
03025         {
03026                 double det;
03027                 register double ted;
03028                 register double *mm=aa->m;
03029                 double a= *(mm++),b= *(mm++);
03030                 double c= *(mm++),d= *(mm);
03031                 det=a*d-b*c;
03032                 if(dtiny(det)) test=0;
03033                 ted=1./det;
03034                 mm=bb->m;
03035                 *(mm++)= d*ted,*(mm++)= -b*ted;
03036                 *(mm++)= -c*ted,*(mm)= a*ted;
03037         }
03038         else if(aa->nr==3)
03039         {
03040                 double det;
03041                 register double ted;
03042                 register double *mm=aa->m;
03043                 double a= *(mm++),b= *(mm++),c= *(mm++);
03044                 double d= *(mm++),e= *(mm++),f= *(mm++);
03045                 double g= *(mm++),h= *(mm++),i= *(mm);
03046                 det=a*e*i-a*h*f-b*d*i+b*g*f+c*d*h-c*g*e;
03047                 if(dtiny(det)) test=0;
03048                 ted=1./det;
03049                 mm=bb->m;
03050 
03051                *(mm++)=(e*i-f*h)*ted,*(mm++)=(c*h-b*i)*ted,*(mm++)=(b*f-e*c)*ted;
03052                 *(mm++)=(f*g-d*i)*ted,*(mm++)=(a*i-g*c)*ted,*(mm++)=(d*c-a*f)*ted;
03053                 *(mm++)=(d*h-g*e)*ted,*(mm++)=(g*b-a*h)*ted,*(mm)=(a*e-d*b)*ted;
03054         }
03055         else
03056         {
03057                 Matrix temp=copy_mx(aa);
03058                 if(gauss_pivot(temp->m,bb->m,aa->nr)==0) test=0;
03059                 close_mx(temp);
03060         }
03061         if(test==0)
03062         {
03063                 e_error("matrix::invert: not invertible, aborting inversion");
03064                 return NULL;
03065         }
03066         return bb;
03067 }
03068 
03069 
03070 
03071 /*-------------------------------------------------------------------------*/
03081 /*--------------------------------------------------------------------------*/
03082 
03083 Matrix
03084 transp_mx(Matrix a)
03085 {
03086         register int nc=a->nc, nr=a->nr;
03087         register double *a0;
03088         register double *b0;
03089         register int i,j;
03090         Matrix b = create_mx(nc,nr);
03091 
03092         if (b == (Matrix)NULL) return b ;
03093         for (i=0,b0=b->m;i<nc;i++)
03094                 for (j=0,a0=a->m+i;j<nr;j++,a0+=nc,b0++)
03095                         *b0 = *a0;
03096         return b;
03097 }
03098 /*-------------------------------------------------------------------------*/
03110 /*--------------------------------------------------------------------------*/
03111 
03112 static int
03113 gauss_pivot(double *ptra, double *ptrc, int n)
03114 /* c(n,n) = a(n,n)^-1 */
03115 {
03116 #define ABS(a) (((a) > 0) ? (a) : -(a))
03117 
03118         register int i,j,k,l;
03119         int maj;
03120         double max,r,t;
03121         double *ptrb;
03122 
03123         ptrb=(double *)calloc(n*n,sizeof(double));
03124         for(i=0;i<n;i++)
03125                 ptrb[i*n+i]= 1.0;
03126 
03127         for (i=1;i <= n;i++)
03128         {
03129                 /* Search max in current column  */
03130                 max = ABS(*(ptra + n*i-n));
03131                 maj = i;
03132                 for (j = i;j <= n;j++)
03133                         if (ABS(*(ptra+n*j+i-n-1)) > max)
03134                         {
03135                                 maj = j;
03136                                 max = ABS(*(ptra+n*j+i-n-1));
03137                         }
03138 
03139                 /* swap lines i and maj */
03140                 if (maj != i)
03141                 {
03142                         for (j = i;j <= n;j++)
03143                         {
03144                                 r = *(ptra+n*maj+j-n-1);
03145                                 *(ptra+n*maj+j-n-1) = *(ptra+n*i+j-n-1);
03146                                 *(ptra+n*i+j-n-1) = r;
03147                         }
03148                         for(l=0;l<n;l++)
03149                         {
03150                                 r = *(ptrb+l*n+maj-1);
03151                                 *(ptrb+l*n+maj-1) = *(ptrb+l*n+i-1);
03152                                 *(ptrb+l*n+i-1) = r;
03153                         }
03154                 }
03155 
03156                /* Subtract line by line */
03157                 for (j = i + 1;j <= n;j++)
03158                 {
03159                         t = (*(ptra+(n+1)*i-n-1));
03160                         if(dtiny(t)) return(0);
03161                         r = (*(ptra+n*j+i-n-1)) / t;
03162                         for(l=0;l<n;l++)
03163                                 *(ptrb+l*n+j-1) -= r * (*(ptrb+l*n+i-1));
03164                         for (k = i;k <= n;k++)
03165                                 *(ptra+n*j+k-n-1) -= r * (*(ptra+n*i+k-n-1));
03166                 }
03167         }
03168 
03169         /* Triangular system resolution */
03170         for(l=0;l<n;l++)
03171                 for (i = n;i >= 1;i--)
03172                 {
03173                         t = (*(ptra+(n+1)*i-n-1));
03174                         if(dtiny(t)) return(0);
03175                         *(ptrc+l+(i-1)*n) = (*(ptrb+l*n+i-1)) / t;
03176                         if (i > 1)
03177                                 for (j = i - 1;j > 0;j--)
03178                                         *(ptrb+l*n+j-1) -= (*(ptra+n*j+i-n-1)) * (*(ptrc+l+(i-1)*n));
03179                 }
03180         free(ptrb);
03181         return(1);
03182 }
03183 /*-------------------------------------------------------------------------*/
03202 /*--------------------------------------------------------------------------*/
03203 
03204 Matrix least_sq_mx(
03205         Matrix  A,
03206         Matrix  B
03207 )
03208 {
03209         Matrix  m1,
03210                         m2,
03211                         m3,
03212                         m4,
03213                         m5 ;
03214 
03215 
03216 
03217         m1 = transp_mx(A) ;
03218         m2 = mul_mx(A, m1) ;
03219         m3 = invert_mx(m2) ;
03220         m4 = mul_mx(B, m1) ;
03221         m5 = mul_mx(m4, m3) ;
03222 
03223         close_mx(m1) ;
03224         close_mx(m2) ;
03225         close_mx(m3) ;
03226         close_mx(m4) ;
03227 
03228         return m5 ;
03229 }
03230 
03231 
03232 /*-------------------------------------------------------------------------*/
03244 /*--------------------------------------------------------------------------*/
03245 
03246 void print_mx(
03247         Matrix  M,
03248         char *  name
03249 )
03250 {
03251         int     i, j ;
03252 
03253         fprintf(stdout, "# matrix %s is [%d x %d]\n", name, M->nr, M->nc) ;
03254         for (j=0 ; j<M->nr ; j++) {
03255                 for (i=0 ; i<M->nc ; i++) {
03256                         fprintf(stdout, "%g\t", M->m[i+j*M->nc]) ;
03257                 }
03258                 fprintf(stdout, "\n") ;
03259         }
03260         fprintf(stdout, "\n") ;
03261 }
03262 
03263 
03264 
03265 /*-------------------------------------------------------------------------*/
03279 /*--------------------------------------------------------------------------*/
03280 
03281 double ipow(double x, int p)
03282 {
03283         double r, recip ;
03284 
03285         /* Get rid of trivial cases */
03286         switch (p) {
03287                 case 0:
03288                 return 1.00 ;
03289 
03290                 case 1:
03291                 return x ;
03292 
03293                 case 2:
03294                 return x*x ;
03295 
03296                 case 3:
03297                 return x*x*x ;
03298 
03299                 case -1:
03300                 return 1.00 / x ;
03301 
03302                 case -2:
03303                 return (1.00 / x) * (1.00 / x) ;
03304         }
03305         if (p>0) {
03306                 r = x ;
03307                 while (--p) r *= x ;
03308         } else {
03309                 r = recip = 1.00 / x ;
03310                 while (++p) r *= recip ;
03311         }
03312         return r;
03313 }
03314 
03315 
03316 /*-------------------------------------------------------------------------*/
03336 /*--------------------------------------------------------------------------*/
03337 
03338 void
03339 ext_free_x(
03340         void    *       pointer,
03341         char    *       filename,
03342         int                     lineno
03343 )
03344 {
03345         malloc_cell     *       curr_cell ;
03346         int                             known_flag = 0 ;
03347 
03348 #ifdef DEBUG_ALLOC
03349         (void)fprintf(stderr,"<- free %p at %s (%d)\n", pointer, filename, lineno);
03350         (void)fflush(stderr) ;
03351 #endif
03352         if (pointer == NULL) {
03353                 e_warning("%s (%d): free requested on NULL pointer\n", filename, lineno);
03354                 return ;
03355         }
03356 
03357         if (malloc_table == NULL) {
03358                 e_warning("%s (%d): free requested on unallocated pointer\n",
03359                                 filename, lineno) ;
03360                 return ;
03361         }
03362 
03363         /*
03364          * Look for the requested pointer to free in our list
03365          */
03366         curr_cell = malloc_table ;
03367         while (curr_cell != NULL) {
03368                 if (curr_cell->pointer == pointer) {
03369                         known_flag = 1 ;
03370                         break ;
03371                 }
03372                 curr_cell = curr_cell->next ;
03373         }
03374         if (known_flag) {
03375                 if (curr_cell->flavour == ALLOC_CELL_SWAP) {
03376                         /*
03377                          * do a swap free: this code is actually a replicate of
03378                          * swap_free(), it saves a function call and a tree search
03379                          * to copy/paste it here.
03380                          */
03381                         total_allocated_SWAP -= curr_cell->size ;
03382                         nalloc_ptrs -- ;
03383                         if (munmap(curr_cell->pointer, curr_cell->size)!=0) {
03384                                 e_error("failed to munmap pointer %s (%d)\n", filename, lineno);
03385                         }
03386                         if (close(curr_cell->swapfd)==-1) {
03387                                 perror("close") ;
03388                                 e_error("closing swap file [%s]", curr_cell->swapfilename);
03389                         }
03390                         if (remove(curr_cell->swapfilename)!=0) {
03391                                 e_error("failed to remove file %s: errno is %d\n",
03392                                                 curr_cell->swapfilename, errno) ;
03393                                 perror("remove") ;
03394                         }
03395                         nswapfiles -- ;
03396                         if ((curr_cell->prev == NULL) && (curr_cell->next == NULL)) {
03397                                 malloc_table = NULL ;
03398                         } else {
03399                                 if (curr_cell->prev == NULL) {
03400                                         /*
03401                                          * First cell becomes malloc_table
03402                                          */
03403                                         malloc_table = curr_cell->next ;
03404                                         malloc_table->prev = NULL ;
03405                                 } else {
03406                                         curr_cell->prev->next = curr_cell->next ;
03407                                 }
03408                                 if (curr_cell->next != NULL) {
03409                                         curr_cell->next->prev = curr_cell->prev ;
03410                                 }
03411                         }
03412                         free(curr_cell) ;
03413                 } else {
03414                         /*
03415                          * do a ram free: this code is actually a replicate of
03416                          * ram_free(), it saves a function call and a tree search
03417                          * to copy/paste it here.
03418                          */
03419                         free(curr_cell->pointer) ;
03420                         total_allocated_RAM -= curr_cell->size ;
03421                         nalloc_ptrs -- ;
03422                         if ((curr_cell->prev == NULL) && (curr_cell->next == NULL)) {
03423                                 malloc_table = NULL ;
03424                         } else {
03425                                 if (curr_cell->prev == NULL) {
03426                                         /*
03427                                          * First cell becomes malloc_table
03428                                          */
03429                                         malloc_table = curr_cell->next ;
03430                                         malloc_table->prev = NULL ;
03431                                 } else {
03432                                         curr_cell->prev->next = curr_cell->next ;
03433                                 }
03434                                 if (curr_cell->next != NULL) {
03435                                         curr_cell->next->prev = curr_cell->prev ;
03436                                 }
03437                         }
03438                         free(curr_cell) ;
03439                 }
03440         } else {
03441                 e_warning("%s (%d) free requested on unallocated pointer\n",
03442                                 filename, lineno) ;
03443         }
03444         return ;
03445 }
03446 
03447 
03448 
03449 /*-------------------------------------------------------------------------*/
03461 /*--------------------------------------------------------------------------*/
03462 
03463 char * fits_getkey(char * line)
03464 {
03465         static char     key[81];
03466         int                             i ;
03467 
03468         if (line==NULL) {
03469 #ifdef DEBUG_FITSHEADER
03470                 printf("fits_getkey: NULL input line\n");
03471 #endif
03472                 return NULL ;
03473         }
03474 
03475         /* Special case: blank keyword */
03476         if (!strncmp(line, "        ", 8)) {
03477                 strcpy(key, "        ");
03478                 return key ;
03479         }
03480         /*
03481          * Sort out special cases: HISTORY, COMMENT and END do not have
03482          * equal signs on the line.
03483          */
03484         if (!strncmp(line, "HISTORY ", 8)) {
03485                 strcpy(key, "HISTORY");
03486                 return key ;
03487         }
03488         if (!strncmp(line, "COMMENT ", 8)) {
03489                 strcpy(key, "COMMENT");
03490                 return key ;
03491         }
03492         if (!strncmp(line, "END ", 4)) {
03493                 strcpy(key, "END");
03494                 return key ;
03495         }
03496 
03497         memset(key, 0, 81);
03498         /* General case: look for the first equal sign */
03499         i=0 ;
03500         while (line[i]!='=' && i<80) i++ ;
03501         if (i>=80) {
03502 #ifdef DEBUG_FITSHEADER
03503                 printf("fits_getkey: cannot find equal sign\n");
03504 #endif
03505                 return NULL ;
03506         }
03507         i-- ;
03508         /* Equal sign found, now backtrack on blanks */
03509         while (line[i]==' ' && i>=0) i-- ;
03510         if (i<=0) {
03511 #ifdef DEBUG_FITSHEADER
03512                 printf("fits_getkey: error backtracking on blanks\n");
03513 #endif
03514                 return NULL ;
03515         }
03516         i++ ;
03517 
03518         /* Copy relevant characters into output buffer */
03519         strncpy(key, line, i) ;
03520         /* Null-terminate the string */
03521         key[i+1] = (char)0;
03522         return key ;
03523 }
03524 
03525 
03526 /*-------------------------------------------------------------------------*/
03539 /*--------------------------------------------------------------------------*/
03540 
03541 char * fits_getcomment(char * line)
03542 {
03543         static char comment[81];
03544         int     i ;
03545         int     from, to ;
03546         int     inq ;
03547 
03548         if (line==NULL) {
03549 #ifdef DEBUG_FITSHEADER
03550                 printf("fits_getcomment: null line in input\n");
03551 #endif
03552                 return NULL ;
03553         }
03554 
03555         /* Special cases: END, HISTORY, COMMENT and blank have no comment */
03556         if (!strncmp(line, "END ", 4)) return NULL ;
03557         if (!strncmp(line, "HISTORY ", 8)) return NULL ;
03558         if (!strncmp(line, "COMMENT ", 8)) return NULL ;
03559         if (!strncmp(line, "        ", 8)) return NULL ;
03560 
03561         memset(comment, 0, 81);
03562         /* Get past the keyword */
03563         i=0 ;
03564         while (line[i]!='=' && i<80) i++ ;
03565         if (i>=80) {
03566 #ifdef DEBUG_FITSHEADER
03567                 printf("fits_getcomment: no equal sign on line\n");
03568 #endif
03569                 return NULL ;
03570         }
03571         i++ ;
03572 
03573         /* Get past the value until the slash */
03574         inq = 0 ;
03575         while (i<80) {
03576                 if (line[i]=='\'')
03577                         inq = !inq ;
03578                 if (line[i]=='/')
03579                         if (!inq)
03580                                 break ;
03581                 i++ ;
03582         }
03583         if (i>=80) {
03584 
03585 
03586 #ifdef DEBUG_FITSHEADER
03587                 printf("fits_getcomment: no slash found on line\n");
03588 #endif
03589                 return NULL ;
03590         }
03591         i++ ;
03592         /* Get past the first blanks */
03593         while (line[i]==' ') i++ ;
03594         from=i ;
03595 
03596         /*
03597          * Now backtrack from the end of the line to the first non-blank
03598          * character.
03599          */
03600         to=79 ;
03601         while (line[to]==' ') to-- ;
03602 
03603         if (to<from) {
03604 #ifdef DEBUG_FITSHEADER
03605                 printf("fits_getcomment: from>to?\n");
03606 #endif
03607                 return NULL ;
03608         }
03609         /* Copy relevant characters into output buffer */
03610         strncpy(comment, line+from, to-from+1);
03611         /* Null-terminate the string */
03612         comment[to-from+1] = (char)0;
03613         return comment ;
03614 }
03615 
03616 
03617 
03618 
03619 /*-------------------------------------------------------------------------*/
03650 /*--------------------------------------------------------------------------*/
03651 
03652 fits_header * fits_read_header(char * filename)
03653 {
03654         mmap_file       *       mm ;
03655         fits_header     *       hdr ;
03656         char                    line[81];
03657         char            *       where ;
03658         char            *       key,
03659                                 *       val,
03660                                 *       com ;
03661 
03662         mm = mmap_open(filename);
03663         if (mm==NULL) return NULL ;
03664 
03665         hdr = fits_header_new();
03666         if (hdr==NULL) {
03667                 mmap_close(mm);
03668                 return NULL ;
03669         }
03670         where = mm->buf ;
03671         while (1) {
03672                 memcpy(line, where, 80);
03673                 line[80] = (char)0;
03674 
03675                 /* Rule out blank lines */
03676                 if (!is_blank_line(line)) {
03677 
03678                         /* Get key, value, comment for current line */
03679                         key = fits_getkey(line);
03680                         val = fits_getvalue(line);
03681                         com = fits_getcomment(line);
03682 
03683                         /* Any NULL response key triggers an error */
03684                         if (key==NULL) {
03685                                 e_error("in header: cannot read key from");
03686                                 e_error("[%s]", line);
03687                                 fits_header_destroy(hdr);
03688                                 hdr = NULL ;
03689                                 break ;
03690                         }
03691 
03692                         /* Append current card to linked-list */
03693                         fits_header_append(hdr, key, val, com, line);
03694                         /* Check for END keyword */
03695                         if (strlen(key)==3)
03696                                 if (key[0]=='E' &&
03697                                         key[1]=='N' &&
03698                                         key[2]=='D')
03699                                         break ;
03700                 }
03701                 /* Move out to next line */
03702                 where += 80 ;
03703                 /* If reaching the end of the file, trigger an error */
03704                 if ((where-mm->buf)>=(mm->size+80)) {
03705                         fits_header_destroy(hdr);
03706                         hdr = NULL ;
03707                         break ;
03708                 }
03709         }
03710         mmap_close(mm);
03711         return hdr ;
03712 }
03713 
03714 
03715 /*-------------------------------------------------------------------------*/
03725 /*--------------------------------------------------------------------------*/
03726 
03727 fits_header * fits_header_new(void)
03728 {
03729         fits_header     *       h ;
03730         h = (fits_header*)llist_create(LISTCOUNT_T_MAX);
03731         return h;
03732 }
03733 
03734 /*-------------------------------------------------------------------------*/
03749 /*--------------------------------------------------------------------------*/
03750 
03751 void fits_header_append(
03752         fits_header     *       hdr,
03753         char    * key,
03754         char    * val,
03755         char    * com,
03756         char    * lin)
03757 {
03758         keytuple        *       k;
03759         llnode_t                *       node ;
03760 
03761         k    = keytuple_new(key, val, com, lin);
03762         node = llnode_create(k);
03763         llist_append((llist_t*)hdr, node);
03764         return ;
03765 }
03766 
03767 
03768 /*---------------------------------------------------------------------------
03769    Function     :       keytuple_new()
03770    In           :       key, value, comment
03771    Out          :       pointer to newly allocated keytuple structure
03772    Job          :       build a new keytuple object
03773    Notice       :       key must be non-NULL, val, com and lin are allowed to be
03774                             NULL or zero-sized character strings.
03775                                 NULL is returned in case of error.
03776  ---------------------------------------------------------------------------*/
03777 
03778 static keytuple * keytuple_new(
03779         char * key,
03780         char * val,
03781         char * com,
03782         char * lin)
03783 {
03784         keytuple        *       k ;
03785 
03786         if (key==NULL) return NULL ;
03787 
03788         /* Allocate space for new structure */
03789         k = malloc(sizeof(keytuple));
03790         /* Hook a copy of the new key */
03791         k->key = strdup2(key) ;
03792         /* Hook a copy of the value if defined */
03793         k->val = NULL ;
03794         if (val!=NULL) {
03795                 if (strlen(val)>0) {
03796                         k->val = strdup2(val);
03797                 }
03798         }
03799         /* Hook a copy of the comment if defined */
03800         k->com = NULL ;
03801         if (com!=NULL) {
03802                 if (strlen(com)>0) {
03803                         k->com = strdup2(com) ;
03804                 }
03805         }
03806         /* Hook a copy of the initial line if defined */
03807         k->lin = NULL ;
03808         if (lin!=NULL) {
03809                 if (strlen(lin)>0) {
03810                         k->lin = strdup2(lin);
03811                 }
03812         }
03813         return k;
03814 }
03815 /*---------------------------------------------------------------------------
03816    Function     :       keytuple_del()
03817    In           :       allocated keytuple object
03818    Out          :       void
03819    Job          :       free the key, val and com fields if they are
03820                             allocated, then free the object itself.
03821    Notice       :       destructor, only to be used with keytuple_new()
03822  ---------------------------------------------------------------------------*/
03823 
03824 static void keytuple_del(keytuple * k)
03825 {
03826         if (k==NULL) return ;
03827         if (k->key) free(k->key);
03828         if (k->val) free(k->val);
03829         if (k->com) free(k->com);
03830         if (k->lin) free(k->lin);
03831         free(k);
03832 }
03833 
03834 
03835 /*---------------------------------------------------------------------------
03836    Function     :       keytuple_cmp()
03837    In           :       two void* pointers, expected to keytuple objects
03838    Out          :       int
03839    Job          :       compare two keytuple objects based on their key field.
03840    Notice       :       returns negative if key1<key2
03841                                                 zero     if key1==key2
03842                                                 positive if key1>key2
03843  ---------------------------------------------------------------------------*/
03844 
03845 static int keytuple_cmp(const void * key1, const void * key2)
03846 {
03847         keytuple * k1,
03848                          * k2 ;
03849 
03850         k1 = (keytuple*)key1 ;
03851         k2 = (keytuple*)key2 ;
03852         return strcmp(k1->key, k2->key);
03853 }
03854 
03855 /*---------------------------------------------------------------------------
03856    Function     :       keytuple_dmp()
03857    In           :       1 keytuple object
03858    Out          :       on stdout, values for key,val,com on one line, value
03859                             of lin on a second line or '--' if NULL.
03860    Job          :       print out a keytuple object on the console
03861    Notice       :       for debugging purposes only
03862  ---------------------------------------------------------------------------*/
03863 
03864 static void keytuple_dmp(keytuple * k)
03865 {
03866         if (!k) return ;
03867         printf("[%s]=[", k->key);
03868         if (k->val) printf("%s", k->val);
03869         printf("]");
03870         if (k->com) printf("/[%s]", k->com);
03871         printf("\n");
03872         return ;
03873 }
03874 
03875 /*-------------------------------------------------------------------------*/
03901 /*--------------------------------------------------------------------------*/
03902 
03903 char * ext_strdup_x(char * s, char * file, int lineno)
03904 {
03905         char * s2 ;
03906 
03907         if (s==NULL) return NULL ;
03908         s2 = ram_malloc(1+(int)strlen(s), file, lineno) ;
03909         strcpy(s2, s) ;
03910         return s2 ;
03911 }
03912 /*-------------------------------------------------------------------------*/
03928 /*--------------------------------------------------------------------------*/
03929 
03930 void destroy_cube_shallow(OneCube * d)
03931 {
03932         int      i ;
03933 
03934     if (d == NULL) return ;
03935 
03936         /* free general plane pointer */
03937         if (d->plane != NULL)
03938                 free(d->plane) ;
03939 
03940     /* free history zone    */
03941         if (d->n_comments>0) {
03942                 for (i=0 ; i<d->n_comments ; i++) {
03943                         free(d->history[i]);
03944                 }
03945                 free(d->history);
03946         }
03947     /* free structure itself    */
03948     free(d) ;
03949 }
03950 
03951 /*-------------------------------------------------------------------------*/
03984 /*--------------------------------------------------------------------------*/
03985 
03986 
03987 poly2d * poly2d_build_from_string(char * s)
03988 {
03989         poly2d  *       p ;
03990         int                     i ;
03991         int                     loop ;
03992         char    *       tok ;
03993         char    *       wstring ;
03994 
03995         if (s==NULL) return NULL ;
03996 
03997         /*
03998          * Count number of values in given string
03999          */
04000         i=0 ;
04001         wstring = strdup2(s);
04002         tok = strtok(wstring, " ");
04003         while (tok!=NULL) {
04004                 i ++ ;
04005                 tok = strtok(NULL, " ");
04006         }
04007         free(wstring);
04008         if (i%3 != 0) {
04009                e_error("in polynomial syntax");
04010                 e_error("the provided string is not made of triplets");
04011                 e_error("[%s]", s);
04012                 return NULL ;
04013         }
04014         i /= 3 ;
04015 
04016         /*
04017          * Allocate space for values
04018          * Parse again the string and assign values into the structure
04019          */
04020 
04021         p = poly2d_allocate(i);
04022         wstring = strdup2(s);
04023         tok = strtok(wstring, " ");
04024         i=0 ;
04025         loop=0 ;
04026         while (tok!=NULL) {
04027                 switch (loop) {
04028                         case 0:
04029                         p->px[i/3] = atoi(tok) ;
04030                         loop++ ;
04031                         break;
04032 
04033                         case 1:
04034                         p->py[i/3] = atoi(tok);
04035                         loop++ ;
04036                         break ;
04037 
04038                         case 2:
04039                         p->c[i/3] = atof(tok);
04040                         loop=0 ;
04041                         break ;
04042                 }
04043                 i++ ;
04044                 tok = strtok(NULL, " ");
04045         }
04046         free(wstring);
04047         return p ;
04048 }
04049 
04050 
04051 /*-------------------------------------------------------------------------*/
04065 /*--------------------------------------------------------------------------*/
04066 
04067 poly2d * poly2d_allocate(int nc)
04068 {
04069         poly2d  *       p ;
04070 
04071         p = malloc(sizeof(poly2d));
04072         p->nc = nc ;
04073         p->px = malloc(nc * sizeof(int));
04074         p->py = malloc(nc * sizeof(int));
04075         p->c  = malloc(nc * sizeof(double));
04076 
04077         return p ;
04078 }
04079 
04080 
04081 
04082 
04083 /*-------------------------------------------------------------------------*/
04093 /*--------------------------------------------------------------------------*/
04094 
04095 void poly2d_free(poly2d * p)
04096 {
04097         if (p==NULL) return ;
04098 
04099         if (p->px != NULL) free(p->px);
04100         if (p->py != NULL) free(p->py);
04101         if (p->c  != NULL) free(p->c);
04102 
04103         free(p);
04104         return ;
04105 }
04106 
04107 /*-------------------------------------------------------------------------*/
04143 /*--------------------------------------------------------------------------*/
04144 
04145 OneImage *
04146 warp_image_generic(
04147     OneImage    *       image_in,
04148         char            *       kernel_type,
04149         poly2d          *       poly_u,
04150         poly2d          *       poly_v
04151 )
04152 {
04153     OneImage    *       image_out ;
04154     int                 i, j, k ;
04155     int                 lx_out, ly_out ;
04156     double              cur ;
04157     double              neighbors[16] ;
04158     double              rsc[8],
04159                                         sumrs ;
04160     double              x, y ;
04161     int                 px, py ;
04162     int                 pos ;
04163     int                 tabx, taby ;
04164     double      *       kernel ;
04165     int                 leaps[16] ;
04166 
04167     if (image_in == NULL) return NULL ;
04168 
04169     /* Generate default interpolation kernel */
04170     kernel = generate_interpolation_kernel(kernel_type) ;
04171     if (kernel == NULL) {
04172         e_error("cannot generate kernel: aborting resampling") ;
04173         return NULL ;
04174     }
04175 
04176     /* Compute new image size   */
04177     lx_out = (int)image_in->lx ;
04178     ly_out = (int)image_in->ly ;
04179 
04180     image_out = new_image(lx_out, ly_out) ;
04181 
04182     /* Pre compute leaps for 16 closest neighbors positions */
04183 
04184     leaps[0] = -1 - image_in->lx ;
04185     leaps[1] =    - image_in->lx ;
04186     leaps[2] =  1 - image_in->lx ;
04187     leaps[3] =  2 - image_in->lx ;
04188 
04189     leaps[4] = -1 ;
04190     leaps[5] =  0 ;
04191     leaps[6] =  1 ;
04192     leaps[7] =  2 ;
04193 
04194     leaps[8] = -1 + image_in->lx ;
04195     leaps[9] =      image_in->lx ;
04196     leaps[10]=  1 + image_in->lx ;
04197     leaps[11]=  2 + image_in->lx ;
04198 
04199     leaps[12]= -1 + 2*image_in->lx ;
04200     leaps[13]=      2*image_in->lx ;
04201     leaps[14]=  1 + 2*image_in->lx ;
04202     leaps[15]=  2 + 2*image_in->lx ;
04203 
04204     /* Double loop on the output image  */
04205     for (j=0 ; j < ly_out ; j++) {
04206         for (i=0 ; i< lx_out ; i++) {
04207             /* Compute the original source for this pixel   */
04208 
04209                         x = poly2d_compute(poly_u, (double)i, (double)j);
04210                         y = poly2d_compute(poly_v, (double)i, (double)j);
04211 
04212             /* Which is the closest integer positioned neighbor?    */
04213             px = (int)x ;
04214                         py = (int)y ;
04215 
04216             if ((px < 1) ||
04217                 (px > (image_in->lx-3)) ||
04218                 (py < 1) ||
04219                 (py > (image_in->ly-3)))
04220                 image_out->data[i+j*lx_out] = (pixelvalue)0.0/0.0 ;
04221             else {
04222                 /* Now feed the positions for the closest 16 neighbors  */
04223                 pos = px + py * image_in->lx ;
04224                 for (k=0 ; k<16 ; k++)
04225                     neighbors[k] =
04226                                         (double)(image_in->data[(int)(pos+leaps[k])]) ;
04227 
04228                 /* Which tabulated value index shall we use?    */
04229                 tabx = (x - (double)px) * (double)(TABSPERPIX) ;
04230                 taby = (y - (double)py) * (double)(TABSPERPIX) ;
04231 
04232                 /* Compute resampling coefficients  */
04233                 /* rsc[0..3] in x, rsc[4..7] in y   */
04234 
04235                 rsc[0] = kernel[TABSPERPIX + tabx] ;
04236                 rsc[1] = kernel[tabx] ;
04237                 rsc[2] = kernel[TABSPERPIX - tabx] ;
04238                 rsc[3] = kernel[2 * TABSPERPIX - tabx] ;
04239                 rsc[4] = kernel[TABSPERPIX + taby] ;
04240                 rsc[5] = kernel[taby] ;
04241                 rsc[6] = kernel[TABSPERPIX - taby] ;
04242                 rsc[7] = kernel[2 * TABSPERPIX - taby] ;
04243 
04244                 sumrs = (rsc[0]+rsc[1]+rsc[2]+rsc[3]) *
04245                         (rsc[4]+rsc[5]+rsc[6]+rsc[7]) ;
04246 
04247                 /* Compute interpolated pixel now   */
04248                 cur =   rsc[4] * (  rsc[0]*neighbors[0] +
04249                                     rsc[1]*neighbors[1] +
04250                                     rsc[2]*neighbors[2] +
04251                                     rsc[3]*neighbors[3] ) +
04252                         rsc[5] * (  rsc[0]*neighbors[4] +
04253                                     rsc[1]*neighbors[5] +
04254                                     rsc[2]*neighbors[6] +
04255                                     rsc[3]*neighbors[7] ) +
04256                         rsc[6] * (  rsc[0]*neighbors[8] +
04257                                     rsc[1]*neighbors[9] +
04258                                     rsc[2]*neighbors[10] +
04259                                     rsc[3]*neighbors[11] ) +
04260                         rsc[7] * (  rsc[0]*neighbors[12] +
04261                                     rsc[1]*neighbors[13] +
04262                                     rsc[2]*neighbors[14] +
04263                                     rsc[3]*neighbors[15] ) ;
04264 
04265                 /* Affect the value to the output image */
04266                 image_out->data[i+j*lx_out] = (pixelvalue)(cur/sumrs) ;
04267                 /* done ! */
04268             }
04269         }
04270     }
04271 
04272     free(kernel) ;
04273     return image_out ;
04274 }
04275 
04276 
04277 /*-------------------------------------------------------------------------*/
04288 /*--------------------------------------------------------------------------*/
04289 
04290 double
04291 sinc(double x)
04292 {
04293     if (fabs(x)<1e-4)
04294         return (double)1.00 ;
04295     else
04296         return ((sin(x * (double)PI_NUMB)) / (x * (double)PI_NUMB)) ;
04297 }
04298 
04299 /*-------------------------------------------------------------------------*/
04323 /*--------------------------------------------------------------------------*/
04324 
04325 double   *
04326 generate_interpolation_kernel(char * kernel_type)
04327 {
04328     double  *   tab ;
04329     int         i ;
04330     double      x ;
04331         double          alpha ;
04332         double          inv_norm ;
04333     int         samples = KERNEL_SAMPLES ;
04334 
04335         if (kernel_type==NULL) {
04336                 tab = generate_interpolation_kernel("tanh") ;
04337         } else if (!strcmp(kernel_type, "default")) {
04338                 tab = generate_interpolation_kernel("tanh") ;
04339         } else if (!strcmp(kernel_type, "sinc")) {
04340                 tab = malloc(samples * sizeof(double)) ;
04341                 tab[0] = 1.0 ;
04342                 tab[samples-1] = 0.0 ;
04343                 for (i=1 ; i<samples ; i++) {
04344                         x = (double)KERNEL_WIDTH * (double)i/(double)(samples-1) ;
04345                         tab[i] = sinc(x) ;
04346                 }
04347         } else if (!strcmp(kernel_type, "sinc2")) {
04348                 tab = malloc(samples * sizeof(double)) ;
04349                 tab[0] = 1.0 ;
04350                 tab[samples-1] = 0.0 ;
04351                 for (i=1 ; i<samples ; i++) {
04352                         x = 2.0 * (double)i/(double)(samples-1) ;
04353                         tab[i] = sinc(x) ;
04354                         tab[i] *= tab[i] ;
04355                 }
04356         } else if (!strcmp(kernel_type, "lanczos")) {
04357                 tab = malloc(samples * sizeof(double)) ;
04358                 for (i=0 ; i<samples ; i++) {
04359                         x = (double)KERNEL_WIDTH * (double)i/(double)(samples-1) ;
04360                         if (fabs(x)<2) {
04361                                 tab[i] = sinc(x) * sinc(x/2) ;
04362                         } else {
04363                                 tab[i] = 0.00 ;
04364                         }
04365                 }
04366         } else if (!strcmp(kernel_type, "hamming")) {
04367                 tab = malloc(samples * sizeof(double)) ;
04368                 alpha = 0.54 ;
04369                 inv_norm  = 1.00 / (double)(samples - 1) ;
04370                 for (i=0 ; i<samples ; i++) {
04371                         x = (double)i ;
04372                         if (i<(samples-1)/2) {
04373                                 tab[i] = alpha + (1-alpha) * cos(2.0*PI_NUMB*x*inv_norm) ;
04374                         } else {
04375                                 tab[i] = 0.0 ;
04376                         }
04377                 }
04378         } else if (!strcmp(kernel_type, "hann")) {
04379                 tab = malloc(samples * sizeof(double)) ;
04380                 alpha = 0.50 ;
04381                 inv_norm  = 1.00 / (double)(samples - 1) ;
04382                 for (i=0 ; i<samples ; i++) {
04383                         x = (double)i ;
04384                         if (i<(samples-1)/2) {
04385                                 tab[i] = alpha + (1-alpha) * cos(2.0*PI_NUMB*x*inv_norm) ;
04386                         } else {
04387                                 tab[i] = 0.0 ;
04388                         }
04389                 }
04390         } else if (!strcmp(kernel_type, "tanh")) {
04391                 tab = generate_tanh_kernel(TANH_STEEPNESS) ;
04392         } else {
04393                 e_error("unrecognized kernel type [%s]: aborting generation",
04394                                 kernel_type) ;
04395                 return NULL ;
04396         }
04397 
04398     return tab ;
04399 }
04400 
04401 
04402 
04403 /*-------------------------------------------------------------------------*/
04415 /*--------------------------------------------------------------------------*/
04416 
04417 #define KERNEL_SW(a,b) tempr=(a);(a)=(b);(b)=tempr
04418 static void reverse_tanh_kernel(double * data, int nn)
04419 {
04420     unsigned long   n,
04421                                         mmax,
04422                                         m,
04423                                         i, j,
04424                                         istep ;
04425     double  wtemp,
04426             wr,
04427             wpr,
04428             wpi,
04429             wi,
04430             theta;
04431     double  tempr,
04432             tempi;
04433 
04434     n = (unsigned long)nn << 1;
04435     j = 1;
04436     for (i=1 ; i<n ; i+=2) {
04437         if (j > i) {
04438             KERNEL_SW(data[j-1],data[i-1]);
04439             KERNEL_SW(data[j],data[i]);
04440         }
04441         m = n >> 1;
04442         while (m>=2 && j>m) {
04443             j -= m;
04444             m >>= 1;
04445         }
04446         j += m;
04447     }
04448     mmax = 2;
04449     while (n > mmax) {
04450         istep = mmax << 1;
04451         theta = 2 * M_PI / mmax;
04452         wtemp = sin(0.5 * theta);
04453         wpr = -2.0 * wtemp * wtemp;
04454         wpi = sin(theta);
04455         wr  = 1.0;
04456         wi  = 0.0;
04457         for (m=1 ; m<mmax ; m+=2) {
04458             for (i=m ; i<=n ; i+=istep) {
04459                 j = i + mmax;
04460                 tempr = wr * data[j-1] - wi * data[j];
04461                 tempi = wr * data[j]   + wi * data[j-1];
04462                 data[j-1] = data[i-1] - tempr;
04463                 data[j]   = data[i]   - tempi;
04464                 data[i-1] += tempr;
04465                 data[i]   += tempi;
04466             }
04467             wr = (wtemp = wr) * wpr - wi * wpi + wr;
04468             wi = wi * wpr + wtemp * wpi + wi;
04469         }
04470         mmax = istep;
04471     }
04472 }
04473 #undef KERNEL_SW
04474 
04475 /*-------------------------------------------------------------------------*/
04496 /*--------------------------------------------------------------------------*/
04497 
04498 #define hk_gen(x,s) (((tanh(s*(x+0.5))+1)/2)*((tanh(s*(-x+0.5))+1)/2))
04499 
04500 double * generate_tanh_kernel(double steep)
04501 {
04502     double  *   kernel ;
04503     double  *   x ;
04504     double      width ;
04505     double      inv_np ;
04506     double      ind ;
04507     int         i ;
04508     int         np ;
04509     int         samples ;
04510 
04511     width   = (double)TABSPERPIX / 2.0 ;
04512     samples = KERNEL_SAMPLES ;
04513     np      = 32768 ; /* Hardcoded: should never be changed */
04514     inv_np  = 1.00 / (double)np ;
04515 
04516     /*
04517      * Generate the kernel expression in Fourier space
04518      * with a correct frequency ordering to allow standard FT
04519      */
04520     x = malloc((2*np+1)*sizeof(double)) ;
04521     for (i=0 ; i<np/2 ; i++) {
04522         ind      = (double)i * 2.0 * width * inv_np ;
04523         x[2*i]   = hk_gen(ind, steep) ;
04524         x[2*i+1] = 0.00 ;
04525     }
04526     for (i=np/2 ; i<np ; i++) {
04527         ind      = (double)(i-np) * 2.0 * width * inv_np ;
04528         x[2*i]   = hk_gen(ind, steep) ;
04529         x[2*i+1] = 0.00 ;
04530     }
04531    /*
04532      * Reverse Fourier to come back to image space
04533      */
04534     reverse_tanh_kernel(x, np) ;
04535 
04536     /*
04537      * Allocate and fill in returned array
04538      */
04539     kernel = malloc(samples * sizeof(double)) ;
04540     for (i=0 ; i<samples ; i++) {
04541         kernel[i] = 2.0 * width * x[2*i] * inv_np ;
04542     }
04543     free(x) ;
04544     return kernel ;
04545 }
04546 
04547 /*-------------------------------------------------------------------------*/
04560 /*--------------------------------------------------------------------------*/
04561 
04562 void show_interpolation_kernel(char * kernel_name)
04563 {
04564         double  *       ker ;
04565         int                     i ;
04566         double          x ;
04567 
04568 
04569         ker = generate_interpolation_kernel(kernel_name) ;
04570         if (ker == NULL)
04571                 return ;
04572 
04573         (void)fprintf(stdout, "# Kernel is %s\n", kernel_name) ;
04574         x = 0.00 ;
04575         for (i=0 ; i<KERNEL_SAMPLES ; i++) {
04576                 (void)fprintf(stdout, "%g %g\n", x, ker[i]) ;
04577                 x += 1.00 / (double)TABSPERPIX ;
04578         }
04579         free(ker) ;
04580         return ;
04581 }
04582 /*-------------------------------------------------------------------------*/
04595 /*--------------------------------------------------------------------------*/
04596 
04597 double
04598 poly2d_compute(
04599         poly2d  *       p,
04600         double          x,
04601         double          y
04602 )
04603 {
04604         double  z ;
04605         int             i ;
04606 
04607         z = 0.00 ;
04608 
04609         for (i=0 ; i<p->nc ; i++) {
04610                 z += p->c[i] * ipow(x, p->px[i]) * ipow(y, p->py[i]) ;
04611         }
04612         return z ;
04613 }
04614 
04615 /*-------------------------------------------------------------------------*/
04625 /*--------------------------------------------------------------------------*/
04626 void fits_header_destroy(fits_header * hdr)
04627 {
04628         llist_t         *       l ;
04629         llnode_t                *       ln ;
04630         keytuple        *       k ;
04631 
04632         l = (llist_t*)hdr ;
04633         for (ln=llist_first(l) ; ln!=0 ; ln=llist_next(l, ln)) {
04634                 k = llnode_get(ln) ;
04635                 keytuple_del(k);
04636         }
04637         llist_destroy_nodes(l);
04638         llist_destroy(l);
04639         return ;
04640 }
04641 /*-------------------------------------------------------------------------*/
04656 /*--------------------------------------------------------------------------*/
04657 char * fits_header_getstr(fits_header * hdr, const char * key)
04658 {
04659         llnode_t                *       ln ;
04660         keytuple        *       k ;
04661 
04662     k = keytuple_new(key, NULL, NULL, NULL) ;
04663     ln = llist_find((llist_t*)hdr, k, keytuple_cmp);
04664     keytuple_del(k);
04665 
04666     if (!ln) return NULL ;
04667     k = llnode_get(ln);
04668     if (!k) return NULL ;
04669         return k->val ;
04670 }
04671 /*-------------------------------------------------------------------------*/
04681 /*--------------------------------------------------------------------------*/
04682 
04683 OneCube *
04684 copy_cube(
04685     OneCube *src_cube
04686 )
04687 {
04688     OneCube     *       dest_cube ;
04689         int                     i ;
04690 
04691     if (src_cube==NULL) return NULL ;
04692     /* allocate data zone, fill up information structure */
04693     dest_cube  = new_cube(  src_cube->lx,
04694                             src_cube->ly,
04695                             src_cube->np) ;
04696     dest_cube ->orig_ptype = src_cube->orig_ptype ;
04697 
04698     /* Then copy data zones from one cube to the other  */
04699     if (src_cube->type==1){
04700         dest_cube->data = (pixelvalue *) calloc(src_cube->np * src_cube->lx * src_cube->ly, sizeof(pixelvalue));
04701         dest_cube->planes= (OneImage *)   calloc(src_cube->np, sizeof(OneImage));
04702         dest_cube->type=1;
04703         for (i=0 ; i<src_cube->np ; i++)
04704         {
04705             fill_new_image(src_cube->lx, src_cube->ly, dest_cube->planes+i, dest_cube->data+(i*src_cube->lx*src_cube->ly));
04706             dest_cube->plane[i] = dest_cube->planes+i;
04707         }
04708         memcpy(dest_cube->data, src_cube->data, (size_t) src_cube->nbpix * sizeof(pixelvalue));
04709     }
04710     else for (i=0 ; i<src_cube->np ; i++) {
04711                 dest_cube->plane[i] = copy_image(src_cube->plane[i]) ;
04712     }
04713     return dest_cube ;
04714 }
04715 
04716 
04717 
04718 /*-------------------------------------------------------------------------*/
04733 /*--------------------------------------------------------------------------*/
04734 
04735 OneCube *
04736 load_cube(const char * filename)
04737 {
04738         OneCube         *       loaded ;
04739 
04740         if (eclipse_is_fits_file(filename)==1) {
04741                 loaded = load_cube_FITS(filename);
04742         } else if (is_ascii_list(filename)==1) {
04743                 loaded = load_cube_ASCII_list(filename);
04744         } else {
04745                 e_error("unsupported format for file [%s]", filename);
04746                 loaded = NULL ;
04747         }
04748         return loaded ;
04749 }
04750 
04751 /*-------------------------------------------------------------------------*/
04762 /*--------------------------------------------------------------------------*/
04763 
04764 int eclipse_is_fits_file(char *filename)
04765 {
04766         FILE    *fp ;
04767         char    *magic ;
04768         int             isfits ;
04769 
04770         if ((fp = fopen(filename, "r"))==NULL) {
04771                 e_error("cannot open file [%s]", filename) ;
04772                 return -1 ;
04773         }
04774 
04775         magic = calloc(FITS_MAGIC_SZ+1, sizeof(char)) ;
04776         (void)fread(magic, 1, FITS_MAGIC_SZ, fp) ;
04777         (void)fclose(fp) ;
04778         magic[FITS_MAGIC_SZ] = (char)0 ;
04779         if (strstr(magic, "SIMPLE")!=NULL)
04780                 isfits = 1 ;
04781         else
04782                 isfits = 0 ;
04783         free(magic) ;
04784         return isfits ;
04785 }
04786 
04787 /*-------------------------------------------------------------------------*/
04796 /*--------------------------------------------------------------------------*/
04797 
04798 OneCube *
04799 load_cube_FITS(char * filename)
04800 {
04801     OneCube     *       loaded_cube ;
04802     OneImage    *       one_plane ;
04803     int                 i ;
04804     cube_info   *       fileinfo ;
04805         mmap_file       *       mm ;
04806         long                    framesize ;
04807         size_t                  pix ;
04808         long                    npix ;
04809         pixelvalue  *   pixbuf ;
04810 
04811     fileinfo = get_cube_info(filename) ;
04812     if (fileinfo == (cube_info*)NULL) {
04813         e_error("in getting FITS information from file %s", filename) ;
04814         return NULL ;
04815     }
04816 
04817     /* Create cube and fill up information fields */
04818     loaded_cube = new_cube(fileinfo->lx, fileinfo->ly, fileinfo->n_im) ;
04819     if (loaded_cube == NULL) {
04820         e_error("cube allocation failure : aborting loading") ;
04821                 free(fileinfo) ;
04822         return NULL ;
04823     }
04824     loaded_cube->orig_ptype = fileinfo->ptype ;
04825 
04826         /* If data formats are compatible, mmap the input file */
04827         if (data_format_compatible(fileinfo->ptype) &&
04828                 (fabs(fileinfo->b_scale - 1.00)<1e-4) &&
04829                 (fabs(fileinfo->b_zero)<1e-4)) {
04830                 mm = mmap_open(filename) ;
04831                 if (mm != NULL) {
04832                         framesize = fileinfo->lx *
04833                                                 fileinfo->ly *
04834                                                 BYTESPERPIXEL(fileinfo->ptype) ;
04835                         for (i=0 ; i<fileinfo->n_im ; i++) {
04836                                 one_plane = malloc(sizeof(OneImage)) ;
04837                                 one_plane->lx = fileinfo->lx ;
04838                                one_plane->ly = fileinfo->ly ;
04839                                 one_plane->nbpix = npix = fileinfo->lx * fileinfo->ly ;
04840                                 one_plane->data = pixbuf =
04841                                         (pixelvalue*)(mm->buf+fileinfo->headersize+i*framesize);
04842 #if BYTE_ORDERING_SCHEME == INTEL
04843                                 for (pix=0 ; pix<npix ; pix++) {
04844                                         swap_bytes(pixbuf+pix, sizeof(pixelvalue));
04845                                 }
04846 #endif
04847                                 /* Register the mmap'ing */
04848                                 one_plane->reg.active = 1 ;
04849                                 one_plane->reg.file_ref = mm->buf ;
04850                                 one_plane->reg.image_ref = pixbuf ;
04851                                 one_plane->reg.mm = mm ;
04852                                 reg_alloc_add(&(one_plane->reg));
04853 
04854                                 loaded_cube->plane[i] = one_plane ;
04855                         }
04856                         free(fileinfo) ;
04857                         return loaded_cube ;
04858                 } else {
04859                         e_warning("mmap failed: loading frame into memory") ;
04860                 }
04861         }
04862 
04863         /* Loop on all planes   */
04864         for (i=0 ; i<fileinfo->n_im ; i++) {
04865                 loaded_cube->plane[i] = load_plane_from_fits(filename, fileinfo, i) ;
04866                 if (loaded_cube->plane[i] == NULL)   {
04867                         e_error("cannot load plane %d: aborting load\n", i+1) ;
04868                         destroy_cube(loaded_cube) ;
04869                         return NULL ;
04870                 }
04871         }
04872 
04873     free(fileinfo) ;
04874     return(loaded_cube) ;
04875 }
04876 
04877 /*-------------------------------------------------------------------------*/
04890 /*--------------------------------------------------------------------------*/
04891 
04892 boolean data_format_compatible(int ptype)
04893 {
04894     int                 compatible ;
04895     float               f ;
04896     double              d ;
04897     unsigned char   *   b ;
04898 
04899     /*
04900      * A word about the pixelvalue type:
04901      * The only restriction put on 'pixelvalue' defined in
04902      * local_types.h, is that the type chosen for that purpose be a
04903      * floating-point type. In ANSI C, this means float, double, or long
04904      * double. The FITS format however only works with 32 bit
04905      * IEEE-compliant floating point and 64 bit IEEE-compliant double
04906      * formats, in Motorola byte ordering (not swapped). The following
04907      * checks will see if indeed a 'pixelvalue' looks like an IEEE float
04908      * in Motorola format or not, or IEEE double.
04909      */
04910 
04911     if ((ptype == BPP_IEEE_FLOAT) &&
04912         (sizeof(pixelvalue)*BITSPERBYTE == 32)) {
04913         /*
04914          * The internal format for pixelvalue is coded on 32 bits
04915          * is it really a float in the IEEE definition?
04916          * Here we choose a floating point number and see if its
04917          * implementation in memory remotely looks like the IEEE
04918          * definition on Motorola machines, namely:
04919          * 2.5e-16 = 0x25901d7d
04920          */
04921         f = (float)2.5e-16 ;
04922         b = (unsigned char*)&f ;
04923 #if BYTE_ORDERING_SCHEME == INTEL
04924                 swap_bytes(b, 4);
04925 #endif
04926 
04927         if ((b[0] == 0x25) &&
04928             (b[1] == 0x90) &&
04929             (b[2] == 0x1d) &&
04930             (b[3] == 0x7d)) {
04931             compatible = 1 ;
04932         } else {
04933             compatible = 0 ;
04934         }
04935    } else if ((ptype == BPP_IEEE_DOUBLE) &&
04936                (sizeof(pixelvalue)*BITSPERBYTE == 64)) {
04937 
04938         /*
04939          * The internal format for pixelvalue is coded on 64 bits
04940          * is it really a double in the IEEE definition?
04941          * Here we choose a double number and see if its implementation
04942          * in memory remotely looks like the IEEE definition on Motorola
04943          * machines, namely:
04944          * 3.14e32 = 0x46aef67970aeecd7
04945          */
04946 
04947         d = (double)3.14e32 ;
04948         b = (unsigned char *)&d ;
04949 #if BYTE_ORDERING_SCHEME == INTEL
04950                 swap_bytes(b, 8);
04951 #endif
04952         if ((b[0] == 0x46) &&
04953             (b[1] == 0xae) &&
04954             (b[2] == 0xf6) &&
04955             (b[3] == 0x79) &&
04956             (b[4] == 0x70) &&
04957             (b[5] == 0xae) &&
04958             (b[6] == 0xec) &&
04959             (b[7] == 0xd7)) {
04960             compatible = 1 ;
04961         } else {
04962             compatible = 0 ;
04963         }
04964     } else {
04965         compatible = 0 ;
04966     }
04967 
04968     return compatible ;
04969 }
04970 
04971 /*-------------------------------------------------------------------------*/
04994 /*--------------------------------------------------------------------------*/
04995 
04996 int
04997 cst_op_image_local(
04998         OneImage        *       image_in,
04999         double                  constant,
05000         int                             operation)
05001 {
05002     int         i ;
05003         double  invlog ;
05004         double  dexp ;
05005         double  invconst ;
05006 
05007     /* Error handling   */
05008         if (image_in == NULL) return -1 ;
05009     if ((fabs((double)constant) < 1e-10) && (operation == '/')){
05010         e_error("division by zero requested in image/constant operation") ;
05011         return -1 ;
05012     }
05013         switch(operation)
05014         {
05015                 case '+':
05016                 for (i=0 ; i<image_in->nbpix ; i++)
05017                         image_in->data[i] =
05018                                 (pixelvalue)((double)image_in->data[i] + constant) ;
05019                 break ;
05020 
05021                 case '-':
05022                 for (i=0 ; i<image_in->nbpix ; i++)
05023                         image_in->data[i] =
05024                                 (pixelvalue)((double)image_in->data[i] - constant) ;
05025                 break ;
05026 
05027                 case '*':
05028                for (i=0 ; i<image_in->nbpix ; i++)
05029                         image_in->data[i] =
05030                                 (pixelvalue)((double)image_in->data[i] * constant) ;
05031                 break ;
05032 
05033                 case '/':
05034                 /* Multiplications are faster than divisions !  */
05035                 invconst = (double)1.0 / constant ;
05036                 for (i=0 ; i<image_in->nbpix ; i++)
05037                         image_in->data[i] =
05038                                 (pixelvalue)((double)image_in->data[i] * invconst) ;
05039                 break ;
05040 
05041                 case 'l':
05042                 invlog = (double)(1.0 / log((double)constant)) ;
05043                 for (i=0 ; i<image_in->nbpix ; i++)
05044                         image_in->data[i] = (pixelvalue)
05045                                                                 log((double)image_in->data[i]*invlog) ;
05046                 break ;
05047 
05048                 case '^':
05049                 dexp = (double)constant ;
05050                 for (i=0 ; i<image_in->nbpix ; i++)
05051                         image_in->data[i] = (pixelvalue)
05052                                                                 pow((double)image_in->data[i], dexp) ;
05053                 break ;
05054 
05055                 case 'e':
05056                 dexp = (double)constant ;
05057                 for (i=0 ; i<image_in->nbpix ; i++)
05058                         image_in->data[i] = (pixelvalue)
05059                                                                 pow(dexp, (double)image_in->data[i]) ;
05060                 break ;
05061 
05062                 default:
05063                 e_error("unrecognized requested operation : aborting") ;
05064                 return -1 ;
05065         }
05066     return 0 ;
05067 }
05068 /*-------------------------------------------------------------------------*/
05091 /*--------------------------------------------------------------------------*/
05092 
05093 OneImage *
05094 shift_image(
05095     OneImage    *       image_in,
05096     double           shift_x,
05097     double           shift_y,
05098     double       *      interp_kernel)
05099 {
05100     OneImage    *       shifted ;
05101     pixelvalue  *       first_pass ;
05102     pixelvalue  *       second_pass ;
05103     int             samples = KERNEL_SAMPLES ;
05104     size_t          i, j ;
05105     double           fx, fy ;
05106     double           rx, ry ;
05107     int             px, py ;
05108     int             tabx, taby ;
05109     double           value ;
05110     size_t          pos ;
05111     register pixelvalue     *   pix ;
05112     register pixelvalue     *   pixint ;
05113     int             mid;
05114     double          norm ;
05115     double       *      ker ;
05116     int                         freeKernel = 1 ;
05117 
05118     /* error handling: test entries */
05119         if (image_in==NULL) return NULL ;
05120 
05121         /* Shifting by a zero offset returns a copy of the input image */
05122         if ((fabs(shift_x)<1e-2) && (fabs(shift_y)<1e-2))
05123                 return copy_image(image_in) ;
05124 
05125         /* See if a kernel needs to be generated */
05126     if (interp_kernel == NULL) {
05127         ker = generate_interpolation_kernel("default") ;
05128         if (ker == NULL) {
05129             e_error("kernel generation failure: aborting resampling") ;
05130             return NULL ;
05131         }
05132     } else {
05133         ker = interp_kernel ;
05134         freeKernel = 0 ;
05135     }
05136 
05137     mid = (int)samples/(int)2 ;
05138     first_pass = calloc(image_in->lx, image_in->ly*sizeof(pixelvalue)) ;
05139     shifted = new_image(image_in->lx, image_in->ly) ;
05140     second_pass = shifted->data ;
05141 
05142     pix = image_in->data ;
05143     for (j=0 ; j<image_in->ly ; j++) {
05144         for (i=1 ; i<image_in->lx-2 ; i++) {
05145             fx = (double)i-shift_x ;
05146             px = (int)fx ;
05147             rx = fx - (double)px ;
05148 
05149             pos = px + j * image_in->lx ;
05150 
05151             if ((px>1) && (px<(image_in->lx-3))) {
05152                 tabx = (int)(fabs((double)mid * rx)) ;
05153                 /*
05154                  * Sum up over 4 closest pixel values,
05155                  * weighted by interpolation kernel values
05156                  */
05157                 value =     (double)pix[pos-1] * ker[mid+tabx] +
05158                             (double)pix[pos] * ker[tabx] +
05159                             (double)pix[pos+1] * ker[mid-tabx] +
05160                             (double)pix[pos+2] * ker[samples-tabx-1] ;
05161                 /*
05162                  * Also sum up interpolation kernel coefficients
05163                  * for further normalization
05164                  */
05165                 norm =      (double)ker[mid+tabx] +
05166                             (double)ker[tabx] +
05167                             (double)ker[mid-tabx] +
05168                             (double)ker[samples-tabx-1] ;
05169                 if (fabs(norm) > 1e-4) {
05170                     value /= norm ;
05171                 }
05172             } else {
05173                 value = 0.0 ;
05174             }
05175             /*
05176              * There may be a problem of rounding here if pixelvalue
05177              * has not enough bits to sustain the accuracy.
05178              */
05179             first_pass[i+j*image_in->lx] = (pixelvalue)value ;
05180         }
05181     }
05182     pixint = first_pass ;
05183     for (i=0 ; i<image_in->lx ; i++) {
05184         for (j=1 ; j<image_in->ly-3 ; j++) {
05185             fy = (double)j - shift_y ;
05186             py = (int)fy ;
05187             ry = fy - (double)py ;
05188             pos = i + py * image_in->lx ;
05189 
05190             taby = (int)(fabs((double)mid * ry)) ;
05191 
05192             if ((py>(int)1) && (py<(image_in->ly-2))) {
05193                 /*
05194                  * Sum up over 4 closest pixel values,
05195                  * weighted by interpolation kernel values
05196                  */
05197                 value = (double)pixint[pos-image_in->lx] * ker[mid+taby] +
05198                         (double)pixint[pos] * ker[taby] +
05199                         (double)pixint[pos+image_in->lx] * ker[mid-taby] +
05200                         (double)pixint[pos+2*image_in->lx]*ker[samples-taby-1];
05201                 /*
05202                  * Also sum up interpolation kernel coefficients
05203                  * for further normalization
05204                  */
05205                 norm =      (double)ker[mid+taby] +
05206                             (double)ker[taby] +
05207                             (double)ker[mid-taby] +
05208                             (double)ker[samples-taby-1] ;
05209 
05210                 if (fabs(norm) > 1e-4) {
05211                     value /= norm ;
05212                 }
05213             } else {
05214                 value = 0.0 ;
05215             }
05216             second_pass[i+j*image_in->lx] = (pixelvalue)value ;
05217         }
05218     }
05219 
05220     free(first_pass) ;
05221     if (freeKernel)
05222         free(ker) ;
05223     return shifted ;
05224 }
05225 
05226 
05227 /*-------------------------------------------------------------------------*/
05239 /*--------------------------------------------------------------------------*/
05240 
05241 int
05242 mul_img_local(
05243         OneImage * im1,
05244         OneImage * im2
05245 )
05246 {
05247         register pixelvalue * p1,
05248                                                 * p2 ;
05249         register int              i ;
05250 
05251         p1 = im1->data ;
05252         p2 = im2->data ;
05253         for (i=0 ; i<im1->nbpix ; i++) {
05254                 *p1++ *= *p2++ ;
05255         }
05256         return 0 ;
05257 }
05258 
05259 
05260 /*-------------------------------------------------------------------------*/
05288 /*--------------------------------------------------------------------------*/
05289 
05290 
05291 #define PIX_SWAP(a,b) { register pixelvalue t=(a);(a)=(b);(b)=t; }
05292 
05293 pixelvalue kth_smallest(pixelvalue a[], int n, int k)
05294 {
05295     register int i,j,l,m ;
05296     register pixelvalue x ;
05297 
05298     l=0 ; m=n-1 ;
05299     while (l<m) {
05300         x=a[k] ;
05301         i=l ;
05302         j=m ;
05303         do {
05304             while (a[i]<x) i++ ;
05305             while (x<a[j]) j-- ;
05306             if (i<=j) {
05307                 PIX_SWAP(a[i],a[j]) ;
05308                 i++ ; j-- ;
05309             }
05310         } while (i<=j) ;
05311         if (j<k) l=i ;
05312         if (k<i) m=j ;
05313     }
05314     return a[k] ;
05315 }
05316 
05317 #undef PIX_SWAP
05318 
05319 
05320 #define PIX_SORT(a,b) { if ((a)>(b)) PIX_SWAP((a),(b)); }
05321 #define PIX_SWAP(a,b) { pixelvalue temp=(a);(a)=(b);(b)=temp; }
05322 #define median_WIRTH(a,n) kth_smallest(a,n,(((n)&1)?((n)/2):(((n)/2)-1)))
05323 
05324 
05325 /*-------------------------------------------------------------------------*/
05339 /*--------------------------------------------------------------------------*/
05340 
05341 pixelvalue
05342 opt_med3(
05343     pixelvalue  *   p
05344 )
05345 {
05346     PIX_SORT(p[0],p[1]) ; PIX_SORT(p[1],p[2]) ; PIX_SORT(p[0],p[1]) ;
05347     return(p[1]) ;
05348 }
05349 
05350 
05351 /*-------------------------------------------------------------------------*/
05365 /*--------------------------------------------------------------------------*/
05366 
05367 pixelvalue
05368 opt_med5(
05369     pixelvalue  *   p
05370 )
05371 {
05372     PIX_SORT(p[0],p[1]) ; PIX_SORT(p[3],p[4]) ; PIX_SORT(p[0],p[3]) ;
05373     PIX_SORT(p[1],p[4]) ; PIX_SORT(p[1],p[2]) ; PIX_SORT(p[2],p[3]) ;
05374     PIX_SORT(p[1],p[2]) ; return(p[2]) ;
05375 }
05376 
05377 
05378 /*-------------------------------------------------------------------------*/
05392 /*--------------------------------------------------------------------------*/
05393 
05394 pixelvalue
05395 opt_med7(
05396     pixelvalue  *   p
05397 )
05398 {
05399     PIX_SORT(p[0], p[5]) ; PIX_SORT(p[0], p[3]) ; PIX_SORT(p[1], p[6]) ;
05400     PIX_SORT(p[2], p[4]) ; PIX_SORT(p[0], p[1]) ; PIX_SORT(p[3], p[5]) ;
05401     PIX_SORT(p[2], p[6]) ; PIX_SORT(p[2], p[3]) ; PIX_SORT(p[3], p[6]) ;
05402     PIX_SORT(p[4], p[5]) ; PIX_SORT(p[1], p[4]) ; PIX_SORT(p[1], p[3]) ;
05403     PIX_SORT(p[3], p[4]) ; return (p[3]) ;
05404 }
05405 /*-------------------------------------------------------------------------*/
05423 /*--------------------------------------------------------------------------*/
05424 
05425 pixelvalue
05426 opt_med9(
05427     pixelvalue  *   p
05428 )
05429 {
05430     PIX_SORT(p[1], p[2]) ; PIX_SORT(p[4], p[5]) ; PIX_SORT(p[7], p[8]) ;
05431     PIX_SORT(p[0], p[1]) ; PIX_SORT(p[3], p[4]) ; PIX_SORT(p[6], p[7]) ;
05432     PIX_SORT(p[1], p[2]) ; PIX_SORT(p[4], p[5]) ; PIX_SORT(p[7], p[8]) ;
05433     PIX_SORT(p[0], p[3]) ; PIX_SORT(p[5], p[8]) ; PIX_SORT(p[4], p[7]) ;
05434     PIX_SORT(p[3], p[6]) ; PIX_SORT(p[1], p[4]) ; PIX_SORT(p[2], p[5]) ;
05435     PIX_SORT(p[4], p[7]) ; PIX_SORT(p[4], p[2]) ; PIX_SORT(p[6], p[4]) ;
05436     PIX_SORT(p[4], p[2]) ; return(p[4]) ;
05437 }
05438 
05439 /*-------------------------------------------------------------------------*/
05457 /*--------------------------------------------------------------------------*/
05458 
05459 
05460 pixelvalue
05461 opt_med25(
05462     pixelvalue  *   p
05463 )
05464 {
05465     PIX_SORT(p[0], p[1]) ;   PIX_SORT(p[3], p[4]) ;   PIX_SORT(p[2], p[4]) ;
05466     PIX_SORT(p[2], p[3]) ;   PIX_SORT(p[6], p[7]) ;   PIX_SORT(p[5], p[7]) ;
05467     PIX_SORT(p[5], p[6]) ;   PIX_SORT(p[9], p[10]) ;  PIX_SORT(p[8], p[10]) ;
05468     PIX_SORT(p[8], p[9]) ;   PIX_SORT(p[12], p[13]) ; PIX_SORT(p[11], p[13]) ;
05469     PIX_SORT(p[11], p[12]) ; PIX_SORT(p[15], p[16]) ; PIX_SORT(p[14], p[16]) ;
05470     PIX_SORT(p[14], p[15]) ; PIX_SORT(p[18], p[19]) ; PIX_SORT(p[17], p[19]) ;
05471     PIX_SORT(p[17], p[18]) ; PIX_SORT(p[21], p[22]) ; PIX_SORT(p[20], p[22]) ;
05472     PIX_SORT(p[20], p[21]) ; PIX_SORT(p[23], p[24]) ; PIX_SORT(p[2], p[5]) ;
05473     PIX_SORT(p[3], p[6]) ;   PIX_SORT(p[0], p[6]) ;   PIX_SORT(p[0], p[3]) ;
05474     PIX_SORT(p[4], p[7]) ;   PIX_SORT(p[1], p[7]) ;   PIX_SORT(p[1], p[4]) ;
05475     PIX_SORT(p[11], p[14]) ; PIX_SORT(p[8], p[14]) ;  PIX_SORT(p[8], p[11]) ;
05476     PIX_SORT(p[12], p[15]) ; PIX_SORT(p[9], p[15]) ;  PIX_SORT(p[9], p[12]) ;
05477     PIX_SORT(p[13], p[16]) ; PIX_SORT(p[10], p[16]) ; PIX_SORT(p[10], p[13]) ;
05478     PIX_SORT(p[20], p[23]) ; PIX_SORT(p[17], p[23]) ; PIX_SORT(p[17], p[20]) ;
05479     PIX_SORT(p[21], p[24]) ; PIX_SORT(p[18], p[24]) ; PIX_SORT(p[18], p[21]) ;
05480     PIX_SORT(p[19], p[22]) ; PIX_SORT(p[8], p[17]) ;  PIX_SORT(p[9], p[18]) ;
05481     PIX_SORT(p[0], p[18]) ;  PIX_SORT(p[0], p[9]) ;   PIX_SORT(p[10], p[19]) ;
05482     PIX_SORT(p[1], p[19]) ;  PIX_SORT(p[1], p[10]) ;  PIX_SORT(p[11], p[20]) ;
05483     PIX_SORT(p[2], p[20]) ;  PIX_SORT(p[2], p[11]) ;  PIX_SORT(p[12], p[21]) ;
05484     PIX_SORT(p[3], p[21]) ;  PIX_SORT(p[3], p[12]) ;  PIX_SORT(p[13], p[22]) ;
05485     PIX_SORT(p[4], p[22]) ;  PIX_SORT(p[4], p[13]) ;  PIX_SORT(p[14], p[23]) ;
05486     PIX_SORT(p[5], p[23]) ;  PIX_SORT(p[5], p[14]) ;  PIX_SORT(p[15], p[24]) ;
05487     PIX_SORT(p[6], p[24]) ;  PIX_SORT(p[6], p[15]) ;  PIX_SORT(p[7], p[16]) ;
05488     PIX_SORT(p[7], p[19]) ;  PIX_SORT(p[13], p[21]) ; PIX_SORT(p[15], p[23]) ;
05489     PIX_SORT(p[7], p[13]) ;  PIX_SORT(p[7], p[15]) ;  PIX_SORT(p[1], p[9]) ;
05490     PIX_SORT(p[3], p[11]) ;  PIX_SORT(p[5], p[17]) ;  PIX_SORT(p[11], p[17]) ;
05491     PIX_SORT(p[9], p[17]) ;  PIX_SORT(p[4], p[10]) ;  PIX_SORT(p[6], p[12]) ;
05492     PIX_SORT(p[7], p[14]) ;  PIX_SORT(p[4], p[6]) ;   PIX_SORT(p[4], p[7]) ;
05493     PIX_SORT(p[12], p[14]) ; PIX_SORT(p[10], p[14]) ; PIX_SORT(p[6], p[7]) ;
05494     PIX_SORT(p[10], p[12]) ; PIX_SORT(p[6], p[10]) ;  PIX_SORT(p[6], p[17]) ;
05495     PIX_SORT(p[12], p[17]) ; PIX_SORT(p[7], p[17]) ;  PIX_SORT(p[7], p[10]) ;
05496     PIX_SORT(p[12], p[18]) ; PIX_SORT(p[7], p[12]) ;  PIX_SORT(p[10], p[18]) ;
05497     PIX_SORT(p[12], p[20]) ; PIX_SORT(p[10], p[20]) ; PIX_SORT(p[10], p[12]) ;
05498 
05499     return (p[12]);
05500 }
05501 
05502 
05503 #undef PIX_SORT
05504 #undef PIX_SWAP
05505 
05506 /*-------------------------------------------------------------------------*/
05521 /*--------------------------------------------------------------------------*/
05522 
05523 pixelvalue median_pixelvalue(pixelvalue * a, int n)
05524 {
05525         pixelvalue      median ;
05526 
05527         switch(n) {
05528                 case 3:
05529                 median = opt_med3(a);
05530                 break ;
05531 
05532                 case 5:
05533                 median = opt_med5(a);
05534                 break ;
05535 
05536                 case 7:
05537                 median = opt_med7(a);
05538                 break ;
05539 
05540                 case 9:
05541                 median = opt_med9(a);
05542                 break ;
05543 
05544                 case 25:
05545                 median = opt_med25(a);
05546                 break ;
05547 
05548                 default:
05549                 median = median_WIRTH(a,n);
05550                 break ;
05551         }
05552         return median;
05553 }
05554 
05555 
05556 
05557 /*-------------------------------------------------------------------------*/
05574 /*--------------------------------------------------------------------------*/
05575 
05576 double image_median_stat(OneImage       *in, double *sigma)
05577 {
05578         double  median_val;
05579         int             i;
05580 
05581         *sigma=0.0;
05582         median_val = (double)get_median_pixelvalue_image(in) ;
05583         for (i=0 ; i<in->nbpix ; i++)
05584                 *sigma+= fabs((double)(in->data[i]-median_val)) ;
05585         *sigma/= (double)in->nbpix ;
05586         return median_val;
05587 }
05588 
05589 /*-------------------------------------------------------------------------*/
05607 /*--------------------------------------------------------------------------*/
05608 
05609 int fillrect_in_image(
05610         OneImage        *       in,
05611         pixelvalue              val,
05612         int                             mini,
05613         int                             maxi,
05614         int                             minj,
05615         int                             maxj)
05616 {
05617     int i;
05618     int j, jl;
05619 
05620         if (in==NULL) return -1 ;
05621         if ((mini<0) ||
05622                 (mini>in->lx) ||
05623                 (maxi>=in->lx) ||
05624                 (minj<0) ||
05625                 (minj>in->ly) ||
05626                 (maxj>=in->ly) ) {
05627                 e_error("incorrect bounds for rectangle fill");
05628                 e_error("got i in [%d-%d] and j in [%d-%d]", mini,maxi,minj,maxj);
05629                 return -1 ;
05630         }
05631     for (j=minj; j<=maxj; j++){
05632         jl = j * in->lx;
05633         for (i=mini; i<=maxi; i++){
05634             in->data[i+jl]= val ;
05635         }
05636     }
05637     return 0;
05638 }
05639 
05640 /*-------------------------------------------------------------------------*/
05651 /*--------------------------------------------------------------------------*/
05652 
05653 OneImage *
05654 image_filter_vertical_median(
05655         OneImage        *       in,
05656         int                             filtsize
05657 )
05658 {
05659         OneImage        *       filt_img=NULL;
05660         int                     col,
05661                                         row;
05662         pixelvalue      *       slit,
05663                                 *       save_column;
05664         int                             maxrow ;
05665         int                             rowdif,
05666                                         f2;
05667 
05668         if ((in==NULL) || (filtsize<1)) return NULL ;
05669 
05670         maxrow = in->ly;
05671         if (maxrow <filtsize){
05672                 e_error("in median filter: filter size is higher than image");
05673                 return NULL;
05674         }
05675 
05676         f2 = filtsize/2;
05677         filt_img = new_image( in->lx, maxrow );
05678         save_column= calloc(in->ly,sizeof(pixelvalue));
05679         slit = calloc(filtsize,sizeof(pixelvalue));
05680         for (col=0;col<in->lx;col++){
05681                 for (row=0;row<in->ly;row++)
05682             save_column[row] = in->data[col+row*in->lx];
05683                 for (row=0;row<in->ly;row++){
05684                         rowdif= f2-row;
05685                         if (rowdif>0)
05686                                 memcpy(slit,save_column,(filtsize-rowdif)*sizeof(pixelvalue));
05687                         else{
05688                                 rowdif= row-in->ly+f2+1;
05689                                 if (rowdif<0)
05690                                         rowdif=0;
05691                                 memcpy(slit,&save_column[row-f2],
05692                                                 (filtsize-rowdif)*sizeof(pixelvalue));
05693                         }
05694                         filt_img->data[col+row*filt_img->lx]=
05695                                 median_pixelvalue(slit,filtsize-rowdif);
05696                 }
05697         }
05698         free(slit);
05699         free(save_column);
05700         return filt_img;
05701 }
05702 
05703 
05704 
05705 /*-------------------------------------------------------------------------*/
05727 /*--------------------------------------------------------------------------*/
05728 
05729 OneCube *
05730 promote_image_to_cube(OneImage * candidate)
05731 {
05732     OneCube * promoted ;
05733 
05734     if (candidate == NULL) return NULL ;
05735     promoted = new_cube(candidate->lx, candidate->ly, 1) ;
05736     promoted->plane[0] = copy_image(candidate) ;
05737     if (promoted->plane[0] == NULL) {
05738         destroy_cube(promoted) ;
05739         return NULL;
05740     }
05741     return promoted ;
05742 }
05743 
05744 
05745 /*-------------------------------------------------------------------------*/
05756 /*--------------------------------------------------------------------------*/
05757 
05758 fits_header * fits_header_default(void)
05759 {
05760         fits_header     *       h ;
05761         h = (fits_header*)llist_create(LISTCOUNT_T_MAX);
05762         fits_header_append(h, "SIMPLE", "T", "Fits format", NULL);
05763         fits_header_append(h, "END", NULL, NULL, NULL);
05764         return h;
05765 }
05766 
05767 /*-------------------------------------------------------------------------*/
05784 /*--------------------------------------------------------------------------*/
05785 
05786 void fits_header_add(
05787         fits_header * hdr,
05788         char    * key,
05789         char    * val,
05790         char    * com,
05791         char    * lin)
05792 {
05793         llist_t         *       l;
05794         keytuple        *       k;
05795         llnode_t                *       node ;
05796         llnode_t                *       end ;
05797 
05798         l = (llist_t*)hdr ;
05799         end  = llist_last(l);
05800         k    = keytuple_new(key, val, com, lin);
05801         node = llnode_create(k);
05802         llist_ins_before(l, node, end);
05803         return ;
05804 }
05805 
05806 /*-------------------------------------------------------------------------*/
05820 /*--------------------------------------------------------------------------*/
05821 void
05822 write_history_in_header(
05823     fits_header  *      fh,
05824     char         **     history,
05825     int                 n_comments
05826 )
05827 {
05828     int i ;
05829 
05830     if (fh==NULL || history==NULL || (n_comments<1)) return ;
05831     for (i=0 ; i<n_comments ; i++) {
05832                 fits_header_add(fh, "HISTORY", history[i], NULL, NULL);
05833     }
05834     return ;
05835 }
05836 /*-------------------------------------------------------------------------*/
05847 /*--------------------------------------------------------------------------*/
05848 int fits_dump_header(
05849         fits_header     *       hdr,
05850         FILE            *       out
05851 )
05852 {
05853         llist_t         *       l ;
05854         llnode_t                *       ln ;
05855         char                    line[81];
05856         int                             n_out ;
05857 
05858         if (hdr==NULL) return -1 ;
05859         if (out==NULL) out=stdout ;
05860 
05861         l = (llist_t*)hdr ;
05862         n_out = 0 ;
05863         for (ln=llist_first(l) ; ln!=0 ; ln=llist_next(l, ln)) {
05864                 /* Make line from information in the node */
05865                 fits_header_makeline(line, ln, 1);
05866                 if ((fwrite(line, 1, 80, out))!=80) {
05867                         fprintf(stderr, "error dumping FITS header");
05868                         return -1 ;
05869                 }
05870                 n_out ++;
05871         }
05872         /* Blank-pad the output */
05873         memset(line, ' ', 80);
05874         while (n_out % 36) {
05875                 fwrite(line, 1, 80, out);
05876                 n_out++ ;
05877         }
05878         return 0 ;
05879 }
05880 
05881 /*-------------------------------------------------------------------------*/
05898 /*--------------------------------------------------------------------------*/
05899 
05900 int append_image_data_to_fits_file(
05901     char            *   filename,
05902     OneImage        *   appended,
05903     int                 pixtype
05904 )
05905 {
05906     FILE    *   out_file ;
05907     BYTE    *   out_buffer ;
05908     int         bytes_per_pixel ;
05909         size_t          to_write,
05910                                 written ;
05911 
05912     /* Error handling : check entries   */
05913     if (appended==NULL || filename==NULL) return -1 ;
05914 
05915     /* Open out file for output in append */
05916         if (!strcmp(filename, "STDOUT")) {
05917                 out_file = stdout ;
05918         } else {
05919                 out_file = fopen(filename, "a") ;
05920         }
05921     if (out_file == NULL) {
05922         e_error("cannot open %s", filename) ;
05923         return -1 ;
05924     }
05925 
05926     /* Allocate a buffer to dump converted data to FITS format  */
05927     bytes_per_pixel = BYTESPERPIXEL(pixtype) ;
05928     out_buffer = malloc(bytes_per_pixel * appended->nbpix) ;
05929     convert_local_to_fits(out_buffer,
05930                           appended->data,
05931                           appended->nbpix,
05932                           pixtype) ;
05933     /* Dump data to disk    */
05934 
05935         to_write = (size_t)((size_t)bytes_per_pixel * appended->nbpix) ;
05936     written = fwrite(out_buffer, sizeof(BYTE), to_write, out_file);
05937         if (to_write != written) {
05938                 e_error("error in writing on disk (disk full?)") ;
05939                e_error("request to write %ld bytes resulted in %ld written to disk",
05940                                 to_write, written) ;
05941         }
05942 
05943         if (out_file != stdout)
05944                 fclose(out_file) ;
05945     free(out_buffer) ;
05946     return 0 ;
05947 }
05948 /*-------------------------------------------------------------------------*/
05963 /*--------------------------------------------------------------------------*/
05964 
05965 void
05966 stuff_with_zeros(
05967     char    *   filename,
05968     size_t              size_before
05969 )
05970 {
05971     FILE    *   stuffed ;
05972     uchar8  *   zeros ;
05973     size_t              added_zeros ;
05974         size_t          written ;
05975 
05976     /* Check silly file sizes   */
05977     if (size_before==0) {
05978         e_error("NULL file size: aborting FITS file zero padding") ;
05979         return ;
05980     }
05981 
05982     /* If the filesize is already a multiple of FITS_BLOCK_SIZE, no job */
05983     if ((size_before % (size_t)FITS_BLOCK_SIZE) == 0)
05984         return ;
05985 
05986     /* How many BYTES should we add to this file ?  */
05987     added_zeros = (size_t)FITS_BLOCK_SIZE -
05988                                   (size_before % (size_t)FITS_BLOCK_SIZE) ;
05989     /* Allocate for a zeroed buffer */
05990     zeros = malloc(added_zeros * sizeof(uchar8)) ;
05991         memset(zeros, 0, added_zeros) ;
05992 
05993     /* Open the file in append  */
05994         if (!strcmp(filename, "STDOUT")) {
05995                 stuffed = stdout ;
05996         } else {
05997                 stuffed = fopen(filename, "a");
05998         }
05999     if (stuffed == NULL) {
06000         e_error("cannot open %s in append", filename) ;
06001         return ;
06002     }
06003 
06004     /* Stuff the file with zeros    */
06005     written = fwrite(zeros, sizeof(uchar8), added_zeros, stuffed) ;
06006    /* OK, quit now */
06007         if (stuffed != stdout)
06008                 fclose(stuffed) ;
06009         if (written != (size_t)added_zeros) {
06010                 e_error("cannot write %ld bytes to file [%s]: aborting (disk full?)",
06011                                 size_before, filename) ;
06012         }
06013     free(zeros) ;
06014 }
06015 
06016 
06017 /*-------------------------------------------------------------------------*/
06027 /*--------------------------------------------------------------------------*/
06028 
06029 int fits_is_boolean(char * s)
06030 {
06031         if (s==NULL) return 0 ;
06032         if (s[0]==0) return 0 ;
06033         if ((int)strlen(s)>1) return 0 ;
06034         if (s[0]=='T' || s[0]=='F') return 1 ;
06035         return 0 ;
06036 }
06037 
06038 
06039 /*-------------------------------------------------------------------------*/
06049 /*--------------------------------------------------------------------------*/
06050 
06051 
06052 int fits_is_int(char * s)
06053 {
06054     regex_t re_int ;
06055     int     status ;
06056 
06057     if (s==NULL) return 0 ;
06058     if (s[0]==0) return 0 ;
06059     if (regcomp(&re_int, &regex_int[0], REG_EXTENDED|REG_NOSUB)!=0) {
06060         e_error(stderr, "internal error: compiling int rule");
06061         exit(-1);
06062     }
06063     status = regexec(&re_int, s, 0, NULL, 0) ;
06064     regfree(&re_int) ;
06065     return (status) ? 0 : 1 ;
06066 
06067 }
06068 
06069 /*-------------------------------------------------------------------------*/
06079 /*--------------------------------------------------------------------------*/
06080 
06081 
06082 int fits_is_float(char * s)
06083 {
06084     regex_t re_float;
06085     int     status ;
06086 
06087     if (s==NULL) return 0 ;
06088     if (s[0]==0) return 0 ;
06089     if (regcomp(&re_float, &regex_float[0], REG_EXTENDED|REG_NOSUB)!=0) {
06090         e_error("internal error: compiling float rule");
06091         exit(-1);
06092     }
06093     status = regexec(&re_float, s, 0, NULL, 0) ;
06094     regfree(&re_float) ;
06095     return (status) ? 0 : 1 ;
06096 }
06097 /*-------------------------------------------------------------------------*/
06107 /*--------------------------------------------------------------------------*/
06108 
06109 int fits_is_complex(char * s)
06110 {
06111     regex_t re_cmp ;
06112     int     status ;
06113 
06114     if (s==NULL) return 0 ;
06115     if (s[0]==0) return 0 ;
06116     if (regcomp(&re_cmp, &regex_cmp[0], REG_EXTENDED|REG_NOSUB)!=0) {
06117         e_error("internal error: compiling complex rule");
06118         exit(-1);
06119     }
06120     status = regexec(&re_cmp, s, 0, NULL, 0) ;
06121     regfree(&re_cmp) ;
06122     return (status) ? 0 : 1 ;
06123 }
06124 
06125 
06126 /*-------------------------------------------------------------------------*/
06141 /*--------------------------------------------------------------------------*/
06142 
06143 void eclipse_keytuple2str(
06144         char * line,
06145         char * key,
06146         char * val,
06147         char * com
06148 )
06149 {
06150         int             i, j ;
06151         int             len ;
06152         int             hierarch = 0 ;
06153         char    cval[81];
06154         char    cval_q[81];
06155         char    ccom[81];
06156 
06157         /* Set the line with zeroes */
06158         memset(line, ' ', 80);
06159         if (key==NULL) return ;
06160 
06161         /* END keyword*/
06162         if (!strcmp(key, "END")) {
06163                 /* Write key and return */
06164                 sprintf(line, "END") ;
06165                 return ;
06166         }
06167         /* HISTORY, COMMENT and blank keywords */
06168         if (!strcmp(key, "HISTORY") ||
06169                 !strcmp(key, "COMMENT") ||
06170                 !strncmp(key, "        ", 8)) {
06171                 /* Write key */
06172                 sprintf(line, "%s ", key);
06173                 if (val==NULL)
06174                         return ;
06175 
06176                 /* There is a value to write, copy it correctly */
06177                 len = strlen(val);
06178                 /*
06179                  * 72 is 80 (FITS line size) minus 8 (sizeof COMMENT or
06180                  * HISTORY)
06181                  */
06182                if (strlen(val)>72) {
06183                         len=72 ;
06184                 }
06185                 strncpy(line+8, val, len);
06186                 return ;
06187         }
06188 
06189         /* Check for NULL values */
06190         if (val==NULL) {
06191                 cval[0]=(char)0;
06192         } else if (strlen(val)<1) {
06193                 cval[0]=(char)0;
06194         } else {
06195                 strcpy(cval, val);
06196         }
06197         /* Check for NULL comments */
06198         if (com==NULL) {
06199                 strcpy(ccom, "no comment");
06200         } else {
06201                 strcpy(ccom, com);
06202         }
06203 
06204         /* Set hierarch flag */
06205         if (!strncmp(key, "HIERARCH", 8)) {
06206                 hierarch ++ ;
06207         }
06208 
06209         /*
06210          * Boolean, int, float or complex
06211          */
06212         if (fits_is_boolean(cval) ||
06213                 fits_is_int(cval) ||
06214                 fits_is_float(cval) ||
06215                 fits_is_complex(cval)) {
06216                 if (hierarch) {
06217                         sprintf(line, "%-29s= %s / %s", key, cval, ccom);
06218                 } else {
06219                         sprintf(line, "%-8.8s= %20s / %-48s", key, cval, ccom);
06220                 }
06221                 return ;
06222         }
06223 
06224         /*
06225      * Blank or NULL values
06226          */
06227 
06228         if (cval[0]==0) {
06229                 if (hierarch) {
06230                         sprintf(line, "%-29s=                    / %s", key, ccom);
06231                 } else {
06232                         sprintf(line, "%-8.8s=                      / %-48s", key, ccom);
06233                 }
06234                 return ;
06235         }
06236         /*
06237          * Can only be a string
06238          * Make simple quotes ['] as double ['']
06239          */
06240 
06241         memset(cval_q, 0, 81);
06242         if (cval[0]=='\'')
06243                 strcpy(cval, fits_pretty_string(cval));
06244         j=0 ;
06245         for (i=0 ; i<strlen(cval) ; i++) {
06246                 if (cval[i]=='\'') {
06247                         cval_q[j]='\'';
06248                         j++ ;
06249                         cval_q[j]='\'';
06250                 } else {
06251                         cval_q[j] = cval[i];
06252                 }
06253                 j++ ;
06254         }
06255         /* sprintf(cval, "'%s'", cval_q); */
06256         if (hierarch) {
06257                 sprintf(line, "%-29s= '%s' / %s", key, cval_q, ccom);
06258         } else {
06259                 sprintf(line, "%-8.8s= '%-18s' / %s", key, cval_q, ccom);
06260         }
06261         return ;
06262 }
06263 
06264 /*-------------------------------------------------------------------------*/
06279 /*--------------------------------------------------------------------------*/
06280 
06281 int fits_header_makeline(
06282         char    * line,
06283         llnode_t * node,
06284         int       conservative)
06285 {
06286         keytuple        *       k ;
06287         char                    blankline[81];
06288         int                             i ;
06289 
06290         k = llnode_get(node);
06291         if (!k) return -1 ;
06292 
06293         /* If a previous line information is there, use it as is */
06294         if (conservative) {
06295                 if (k->lin != NULL) {
06296                         memcpy(line, k->lin, 80);
06297                         line[80]=(char)0;
06298                         return 0 ;
06299                 }
06300         }
06301         /* Got to build keyword from scratch */
06302         memset(blankline, 0, 81);
06303         eclipse_keytuple2str(blankline, k->key, k->val, k->com);
06304         memset(line, ' ', 80);
06305         line[80]=(char)0;
06306         for (i=0 ; i<strlen(blankline) ; i++) {
06307                 line[i] = blankline[i];
06308         }
06309         return 0;
06310 }
06311 /*-------------------------------------------------------------------------*/
06334 /*--------------------------------------------------------------------------*/
06335 
06336 void
06337 convert_local_to_fits(
06338         BYTE                    *p_dest,
06339         pixelvalue              *p_source,
06340         ulong32                             nbpix,
06341         int                             pixel_type )
06342 {
06343         int                     i ;
06344         short16         shortintpix ;
06345         long32          longintpix ;
06346         float           floatpix ;
06347         double          doublepix ;
06348         BYTE            *mem_dest ;
06349 
06350         /* Error handling: test entries */
06351 
06352         if (p_dest==NULL || p_source==NULL) return ;
06353         if ((nbpix<1)||(nbpix > MAX_COLUMN_NUMBER * MAX_LINE_NUMBER)) return ;
06354 
06355         mem_dest = p_dest ;
06356     switch(pixel_type) {
06357 
06358         case BPP_8_UNSIGNED:
06359                 /* No need to care about endian : we are dealing with bytes     */
06360                 for (i=0 ; i<nbpix ; i++) {
06361                         if (p_source[i] > (pixelvalue)UCHAR_MAX) {
06362                                 /* Clip the biggest number to UCHAR_MAX */
06363                                 p_dest[i] = (BYTE)UCHAR_MAX ;
06364                         } else {
06365                                 if (p_source[i] < (pixelvalue)0) {
06366                                         /* Clip the smallest to 0   */
06367                                         mem_dest[i] = (BYTE)0 ;
06368                                 } else {
06369                                         /* No clipping : just cast to smallest neighbor  */
06370                                         mem_dest[i] = (BYTE)floor((double)p_source[i]+0.5) ;
06371                                 }
06372                         }
06373                 }
06374                 break ;
06375 
06376         case BPP_16_SIGNED:
06377                 for (i=0 ; i<nbpix ; i++) {
06378                         if (p_source[i] > (pixelvalue)SHRT_MAX) {
06379                                 /* Clip the biggest number to SHRT_MAX */
06380                                 shortintpix = (short16)(SHRT_MAX);
06381                         } else {
06382                                 if (p_source[i] < (pixelvalue)SHRT_MIN) {
06383                                         /* Clip the smallest to SHRT_MIN   */
06384                                         shortintpix = (short16)SHRT_MIN ;
06385                                 } else {
06386                                         /* No clipping : just cast to closest neighbor */
06387                                         shortintpix = (short16)floor((double)p_source[i]+0.5) ;
06388                                 }
06389                         }
06390 
06391 #if (BYTE_ORDERING_SCHEME == INTEL)
06392                         shortintpix = swap_2_bytes(shortintpix) ;
06393 #endif
06394                         memcpy(mem_dest, &shortintpix, 2) ;
06395                         mem_dest += 2 ;
06396                 }
06397         break ;
06398 
06399         case BPP_32_SIGNED:
06400                 for (i=0 ; i<nbpix ; i++) {
06401                         if (p_source[i] > (pixelvalue)LONG_MAX) {
06402                                 /* Clip the biggest number to LONG_MAX */
06403                                 longintpix = (long32)LONG_MAX ;
06404                         } else {
06405                                 if (p_source[i] < (pixelvalue)LONG_MIN) {
06406                                         /* Clip the smallest to LONG_MIN   */
06407                                         longintpix = (long32)LONG_MIN ;
06408                                 } else {
06409                                         /* No clipping : just cast  */
06410                                         longintpix = (long32)floor((double)p_source[i]+0.5) ;
06411                                 }
06412                         }
06413 
06414 #if (BYTE_ORDERING_SCHEME == INTEL)
06415                         longintpix = swap_4_bytes(longintpix) ;
06416 #endif
06417 
06418                         memcpy(mem_dest, &longintpix, 4) ;
06419                         mem_dest += 4 ;
06420                 }
06421         break ;
06422 
06423         case BPP_IEEE_FLOAT:
06424                 for (i=0 ; i<nbpix ; i++) {
06425                         /*
06426                          * This part is tricky: how to convert a pixelvalue which is
06427                          * a float or a double, into 4 swappable bytes?
06428                          * The algorithm is the following:
06429                          * 1. convert the pixel to a float (cast operator may loose
06430                          *    significant bits).
06431                          * 2. transfer the bits of the float to a long
06432                          * 3. swap the long if needed
06433                          * 4. write the bits to an output buffer.
06434                          */
06435                         floatpix = (float)p_source[i] ;
06436                         longintpix = *((long*)&floatpix) ;
06437 #if (BYTE_ORDERING_SCHEME == INTEL)
06438                 /* 32 bit data bytes swapping   */
06439                         longintpix = swap_4_bytes(longintpix) ;
06440 #endif
06441                         memcpy(mem_dest, &longintpix, 4) ;
06442                         mem_dest += 4 ;
06443                 }
06444         break ;
06445 
06446 
06447                 case BPP_IEEE_DOUBLE:
06448                 for (i=0 ; i<nbpix ; i++) {
06449                         doublepix = (double)p_source[i] ;
06450 #if (BYTE_ORDERING_SCHEME == INTEL)
06451                         /* swap_8_bytes((BYTE*)(&doublepix)) ; */
06452                         swap_bytes(&doublepix, 8) ;
06453 #endif
06454                         memcpy(mem_dest, &doublepix, 8) ;
06455                         mem_dest += 8 ;
06456                 }
06457         break ;
06458 
06459 
06460         default:
06461         e_error("unrecognized pixel type in output: aborting conversion") ;
06462         break ;
06463     }
06464 }
06465 
06466 
06467 /*-------------------------------------------------------------------------*/
06486 /*--------------------------------------------------------------------------*/
06487 
06488 void
06489 save_cube_to_fits(
06490     OneCube     *       to_save,
06491     char        *       filename
06492 )
06493 {
06494         fits_header     *       fh ;
06495         FILE            *       out ;
06496         char                    cval[80];
06497     int                         i ;
06498     size_t                      total_size ;
06499 
06500         /* Sanity checks */
06501     if ((to_save == NULL) || (filename==NULL)) return ;
06502     if ((to_save->lx > MAX_COLUMN_NUMBER) ||
06503         (to_save->lx < 1) ||
06504         (to_save->ly > MAX_LINE_NUMBER) ||
06505         (to_save->ly < 1) ||
06506         (to_save->np > MAX_IMAGE_NUMBER) ||
06507         (to_save->np < 1) ||
06508         (to_save->nbpix != ((size_t)to_save->lx *
06509                             (size_t)to_save->ly *
06510                             (size_t)to_save->np))) {
06511         e_error("error in entry cube : inconsistent values !") ;
06512                 e_error("size is [%ld x %ld x %ld], nbpix is %ld",
06513                                 to_save->lx,
06514                                 to_save->ly,
06515                                 to_save->np,
06516                                 to_save->nbpix) ;
06517         return ;
06518     }
06519 
06520         /* Make a default FITS header for this cube */
06521         fh = fits_header_default();
06522         /* BITPIX */
06523         sprintf(cval, "%d", to_save->orig_ptype);
06524         fits_header_add(fh, "BITPIX", cval, "Bits per pixel", NULL);
06525 
06526         /* NAXIS */
06527         if (to_save->np > 1) {
06528                 fits_header_add(fh, "NAXIS", "3", "File dimension", NULL);
06529         } else {
06530                 fits_header_add(fh, "NAXIS", "2", "File dimension", NULL);
06531         }
06532         /* NAXIS1 NAXIS2 (NAXIS3) */
06533         sprintf(cval, "%d", to_save->lx);
06534         fits_header_add(fh, "NAXIS1", cval, "Size in x", NULL);
06535         sprintf(cval, "%d", to_save->ly);
06536         fits_header_add(fh, "NAXIS2", cval, "Size in y", NULL);
06537         if (to_save->np>1) {
06538                 sprintf(cval, "%d", to_save->np);
06539                 fits_header_add(fh, "NAXIS3", cval, "Number of planes", NULL);
06540         }
06541     /*
06542      * Adding BSCALE and BZERO keywords for Workcid compatibility
06543      * eclipse works internally without scaling and offset, so the
06544      * values shall always be BSCALE=1 and BZERO=0.
06545      */
06546         fits_header_add(fh, "BSCALE", "1.0", "pixel scale factor", NULL);
06547         fits_header_add(fh, "BZERO",  "0.0", "pixel offset", NULL);
06548         /* Add eclipse signature */
06549         fits_header_add(fh, "ECLIPSE", "1", "created by eclipse", NULL);
06550         fits_header_add(fh, "ORIGIN", "eclipse", "created by eclipse", NULL);
06551     /* Write the eclipse history string into header */
06552     write_history_in_header(fh, to_save->history, to_save->n_comments) ;
06553 
06554     /* Ouput header to file */
06555         if (!strcmp(filename, "STDOUT")) {
06556                 out = stdout ;
06557         } else {
06558                 out = fopen(filename, "w");
06559         }
06560         fits_dump_header(fh, out);
06561         if (out!=stdout)
06562                 fclose(out);
06563         fits_header_destroy(fh);
06564 
06565    /* Convert planes one by one and copy them into Fits structure  */
06566     for (i=0 ; i<to_save->np ; i++) {
06567                 if (to_save->np>1)
06568                         compute_status("converting plane", i, to_save->np, 3) ;
06569         if (append_image_data_to_fits_file( filename,
06570                                                                                         to_save->plane[i],
06571                                                                                         to_save->orig_ptype)!=0){
06572             e_error("cannot append plane %d to file [%s]: aborting save",
06573                                         i+1, filename) ;
06574             return ;
06575         }
06576     }
06577 
06578     /*
06579          * Now stuff the rest of the file until a multiple of FITS_BLOCK_SIZE
06580      * is reached (see FITS format spec).
06581      */
06582 
06583     total_size = BYTESPERPIXEL(to_save->orig_ptype) ;
06584     total_size *= to_save->nbpix ;
06585     stuff_with_zeros(filename, total_size) ;
06586 
06587 }
06588 
06589 /*-------------------------------------------------------------------------*/
06619 /*--------------------------------------------------------------------------*/
06620 
06621 void
06622 save_image_to_fits(
06623     OneImage    *       to_save,
06624     char        *       filename,
06625     int                 pixel_type
06626 )
06627 {
06628     OneCube *   saved_cube ;
06629 
06630     if ((to_save==NULL) || (filename==NULL)) return ;
06631 
06632     saved_cube = promote_image_to_cube(to_save) ;
06633     if (saved_cube == NULL) {
06634         e_error("cannot promote image to cube: aborting save image") ;
06635         return ;
06636     }
06637     saved_cube->orig_ptype = pixel_type ;
06638     save_cube_to_fits(saved_cube, filename) ;
06639     destroy_cube(saved_cube) ;
06640 }
06641 
06642 
06643 
06644 
06645 /*---------------------------------------------------------------------------*/
06659 /*---------------------------------------------------------------------------*/
06660 static int threshold_one_dim(
06661                 OneImage                *       in,
06662                 pixelvalue                      threshold,
06663                 pixelvalue              *       line,
06664                 orientation_t           orientation,
06665                 pixelvalue                      val)
06666 {
06667         int                     i,
06668                                         j ;
06669         int                             err = 0 ;
06670 
06671         /**************/
06672         /* DEBUG CASE */
06673         if (debug_active()>=ARC_DEBUG_LEVEL)
06674                 e_comment(0, "threshold_one_dim: thr=%f", (float)threshold) ;
06675         /**************/
06676 
06677         threshold=threshold+1;
06678         switch(orientation) {
06679                 case VERTICAL:
06680                 for (i=0 ; i<in->lx ; i++)
06681                         if (line[i] < threshold) {
06682                                 for (j=0 ; j<in->ly ; j++) in->data[j*in->lx+i] = val ;
06683                         }
06684                 break ;
06685 
06686                 case HORIZONTAL:
06687                 for (j=0 ; j<in->ly ; j++)
06688                         if (line[j] < threshold) {
06689                                 for (i=0 ; i<in->lx ; i++) in->data[j*in->lx+i] = val ;
06690                         }
06691                 break;
06692 
06693                 default:
06694                 e_error("threshold_one_dim: unknown orientation %d",
06695                                 (int)orientation) ;
06696                 err = -1 ;
06697         }
06698 
06699         /***************/
06700         /* DEBUG CASE  */
06701 
06702        if ((debug_active()>=ARC_DEBUG_LEVEL) && !err)
06703                 save_image_to_fits(in, "sub1d.fits", ARC_PARAM_PTYPE) ;
06704         /***************/
06705 
06706         return err ;
06707 }
06708 
06709 /*-------------------------------------------------------------------------*/
06723 /*--------------------------------------------------------------------------*/
06724 
06725 pixel_map *
06726 new_pixelmap(
06727     int     lx,
06728     int     ly
06729 )
06730 {
06731     pixel_map * p ;
06732 
06733     if ((lx<=0)||(lx>MAX_COLUMN_NUMBER)||(ly<=0)||(ly>MAX_LINE_NUMBER)) {
06734         e_error("size out of bounds: cannot create pixel map") ;
06735         e_error("requested size is [%d x %d]", lx, ly) ;
06736         return NULL ;
06737     }
06738     p = malloc(sizeof(pixel_map)) ;
06739     p->lx = lx ;
06740     p->ly = ly ;
06741     p->nbpix = p->ngoodpix = lx * ly ;
06742     p->data = malloc(lx * ly * sizeof(binpix));
06743     /* Set all pixel values to 1 */
06744     memset(p->data, PIXELMAP_1, p->nbpix);
06745     return p ;
06746 }
06747 
06748 
06749 /*-------------------------------------------------------------------------*/
06765 /*--------------------------------------------------------------------------*/
06766 
06767 pixel_map *
06768 thresh_image_to_pixelmap(
06769     OneImage    *       in,
06770     double              lo_cut,
06771     double              hi_cut
06772 )
06773 {
06774     pixel_map    *      p ;
06775     int                 i ;
06776 
06777         if (in==NULL) return NULL ;
06778     p = new_pixelmap(in->lx, in->ly) ;
06779     /* Loop on all pixels   */
06780     for (i=0 ; i<p->nbpix ; i++) {
06781         if (in->data[i] < lo_cut) {
06782             p->data[i] = PIXELMAP_0 ;
06783             p->ngoodpix -- ;
06784         } else if (in->data[i] > hi_cut) {
06785             p->data[i] = PIXELMAP_0 ;
06786             p->ngoodpix -- ;
06787         }
06788     }
06789     return p ;
06790 }
06791 
06792 /*-------------------------------------------------------------------------*/
06802 /*--------------------------------------------------------------------------*/
06803 
06804 pixelvalue
06805 get_median_pixelvalue_image(OneImage * in)
06806 {
06807         OneImage        *       clone ;
06808         pixelvalue              median ;
06809 
06810         clone = copy_image(in) ;
06811         if (clone == NULL) {
06812                 return 0.00 ;
06813         }
06814         median = median_pixelvalue(clone->data, clone->nbpix) ;
06815         destroy_image(clone) ;
06816         return median ;
06817 }
06818 
06819 
06820 /*-------------------------------------------------------------------------*/
06831 /*--------------------------------------------------------------------------*/
06832 
06833 OneImage *
06834 image_filter_horizontal_median(
06835         OneImage        *       in,
06836         int                             filtsize )
06837 {
06838         OneImage        *       filt_img=NULL;
06839         int                     col,
06840                                         row;
06841         pixelvalue      *       save_slit;
06842         int                             maxcol;
06843         int                             f2 ;
06844         int                             coldif;
06845 
06846         maxcol = in->lx;
06847         f2 = filtsize/2;
06848         if (maxcol <filtsize){
06849                 e_error("in median filter: filter is larger than image");
06850                 return NULL;
06851         }
06852         filt_img =  new_image( maxcol,in->ly );
06853         save_slit = calloc(filtsize,sizeof(pixelvalue));
06854         for (row=0;row<in->ly;row++){
06855                 for (col=0;col<maxcol;col++){
06856                         coldif= f2-col;
06857                         if (coldif>0)
06858                                 memcpy(save_slit,&(in->data[0+row*in->lx]),
06859                                                 (filtsize-coldif)*sizeof(pixelvalue));
06860                         else{
06861                                 coldif= col-in->lx+f2+1;
06862                                 if (coldif<0)
06863                                         coldif=0;
06864                                         memcpy(save_slit,&(in->data[col-f2+row*in->lx]),
06865                                                 (filtsize-coldif)*sizeof(pixelvalue));
06866                         }
06867                         filt_img->data[col+row*filt_img->lx]=
06868                                 median_pixelvalue(save_slit,filtsize-coldif);
06869                 }
06870         }
06871         free(save_slit);
06872         return filt_img;
06873 }
06874 
06875 
06876 
06877 /*-------------------------------------------------------------------------*/
06895 /*--------------------------------------------------------------------------*/
06896 
06897 
06898 
06899 OneImage *
06900 collapse_image_median(
06901         OneImage        *       in,
06902         int                             direction,
06903         int                             discard_lo,
06904         int                             discard_hi
06905 )
06906 {
06907         OneImage        *       collapsed ;
06908         pixelvalue      *       line ;
06909         int                             i, j ;
06910         int                             width ;
06911 
06912         if (in==NULL) return NULL ;
06913         if (direction == 1) {
06914                 /* Collapsing the image in the x direction */
06915                 if ((discard_lo+discard_hi)>=in->ly) {
06916                         e_error("discard bounds: %d+%d >= %d",
06917                                         discard_lo, discard_hi, in->ly);
06918                         return NULL ;
06919                 }
06920                 collapsed = new_image(1, in->ly);
06921                 width = in->lx - discard_lo - discard_hi ;
06922                 line = malloc(width * sizeof(pixelvalue)) ;
06923                 for (j=0 ; j<in->ly ; j++) {
06924                         for (i=0 ; i<width ; i++) {
06925                                 line[i] = in->data[(i+discard_lo)+j*in->lx];
06926                         }
06927                         collapsed->data[j] = median_pixelvalue(line, width);
06928                 }
06929                 free(line);
06930         } else if (direction==0) {
06931                 /* Collapsing the image in the y direction */
06932                 if ((discard_lo+discard_hi)>=in->lx) {
06933                         e_error("discard bounds: %d+%d >= %d",
06934                                         discard_lo, discard_hi, in->lx);
06935                         return NULL ;
06936                 }
06937                 collapsed = new_image(in->lx, 1);
06938                 width = in->ly - discard_lo - discard_hi ;
06939                 line = malloc(width * sizeof(pixelvalue)) ;
06940                 for (i=0 ; i<in->lx ; i++) {
06941                         for (j=0 ; j<width ; j++) {
06942                                 line[j] = in->data[i+(j+discard_lo)*in->lx] ;
06943                         }
06944                         collapsed->data[i] = median_pixelvalue(line, width);
06945                 }
06946                 free(line);
06947         } else {
06948                 e_error("unknown direction for collapsing: [%d]", direction);
06949                 collapsed = NULL ;
06950         }
06951         return collapsed ;
06952 }
06953 
06954 
06955 /*-------------------------------------------------------------------------*/
06965 /*--------------------------------------------------------------------------*/
06966 
06967 void destroy_pixelmap(pixel_map * p)
06968 {
06969     if (p == NULL) return ;
06970     if (p->data != NULL) {
06971         free(p->data) ;
06972     }
06973     free(p) ;
06974     return ;
06975 }
06976 
06977 /*--------------------------------------------------------------------------*/
06986 /*--------------------------------------------------------------------------*/
06987 label_image_t * new_label_image(
06988                 size_t  lx,
06989                 size_t  ly)
06990 {
06991         label_image_t   *       lab ;
06992         size_t                  nbpix ;
06993 
06994         /* Check the requested size */
06995         if ((lx <= 0) || (lx > MAX_COLUMN_NUMBER) ||
06996                 (ly <= 0) || (ly > MAX_LINE_NUMBER)) {
06997         e_error("size out of bounds: cannot create label image") ;
06998         e_error("requested size is [%d x %d]", lx, ly) ;
06999         return NULL ;
07000         }
07001 
07002         /* Set up nbpix */
07003         nbpix = lx * ly ;
07004 
07005         /* Memory allocation */
07006         lab = (label_image_t*)calloc(1, sizeof(label_image_t)) ;
07007 
07008         /* Fill the fields */
07009         lab->lx = lx ;
07010         lab->ly = ly ;
07011         lab->nrealobj = 0 ;
07012         lab->nobjtot = 0 ;
07013 
07014         /* Pixels allocation */
07015     lab->buf = malloc(nbpix * sizeof(label_pix)) ;
07016     if (lab == NULL) {
07017         e_error("memory allocation failure: aborting label image creation") ;
07018         free(lab) ;
07019         return NULL ;
07020     }
07021 
07022         /* Return */
07023     return lab ;
07024 }
07025 /*--------------------------------------------------------------------------*/
07041 /*--------------------------------------------------------------------------*/
07042 static label_pix search_n_classify(
07043                 label_image_t   *       l,
07044                 label_pix                       first_class,
07045                 cleq_tab_t              *       cl)
07046 {
07047         /* row and col   */
07048         size_t                  i,
07049                                         j ;
07050         /* row*col combinations for faster indexing in image */
07051         size_t                  jl ;
07052         /* flag signalling class equiavlence table not saturated */
07053         int                             table_ok ;
07054         /* intermediate class counter */
07055         label_pix       curclass ;
07056         /* pointer to destinatation pixel */
07057         label_pix       *       pdest ;
07058         /* pixel neighbourhood definition */
07059         label_pix               p0_1,
07060                                 p_10,
07061                                         p10,
07062                                         p01 ;
07063         /* labels assigned within neighbourhood */
07064         label_pix               lb0_1,
07065                                 lb_10,
07066                                         lb00,
07067                                         lb10,
07068                                         lb01 ;
07069 
07070         /* Initialise */
07071         curclass = first_class ;
07072         table_ok = TRUE ;
07073 
07074 
07075         for (j=0 ; j<l->ly ; j++) {
07076                 jl = j * l->lx ;
07077         for (i=0 ; i<l->lx ; i++) {
07078                 pdest = &l->buf[i+jl] ;
07079                         if (*pdest == CNCTCMP_BLACK) continue ;
07080                 if (*pdest != CNCTCMP_WHITE) {
07081                                 if (curclass <= *pdest) curclass = *pdest + 1 ;
07082                                 if (table_ok) table_ok = add_equiv_cl(cl, *pdest, *pdest) ;
07083                                 if (table_ok) table_ok = add_equiv_cl(cl, *pdest, *pdest) ;
07084                                 continue ;
07085                         }
07086                         get4nbs(l, i, j, &p0_1, &p_10, &p10, &p01) ;
07087                         lb_10 = lb0_1 = lb10 = lb01 = 0 ;
07088 
07089                         /* check defined neighbours*/
07090                 switch(p_10) {
07091                         case CNCTCMP_BLACK:
07092                                 /* skip */
07093                         break ;
07094 
07095                         case CNCTCMP_WHITE:
07096                                 /* preceding point should have been labelled */
07097                 e_error("labeling order error at p_10(%d,%d)=%d", i, j, p_10) ;
07098                         return 0 ;
07099 
07100                                 default:
07101                                 lb_10 = p_10 ;
07102                 }
07103 
07104                 switch(p0_1) {
07105                         case CNCTCMP_BLACK:
07106                                 /*skip*/
07107                         break ;
07108 
07109                                 case CNCTCMP_WHITE:
07110                                 /*preceding point should have been labelled*/
07111                         e_error("labeling order error at p0_1(%d,%d)", i, j) ;
07112                         return 0 ;
07113 
07114                                 default:
07115                                 if (p0_1 != lb_10) lb0_1 = p0_1 ;
07116                 }
07117 
07118                 switch(p01) {
07119                                 /*unlabelled*/
07120                         case CNCTCMP_BLACK:
07121                         case CNCTCMP_WHITE:
07122                                 break ;
07123                         default:
07124                                 if ((p01 != lb_10) && (p01 !=lb0_1)) lb01 = p01 ;
07125                 }
07126 
07127                 switch(p10){
07128                                 /*unlabelled*/
07129                         case CNCTCMP_BLACK:
07130                         case CNCTCMP_WHITE:
07131                                 break ;
07132                         default:
07133                                 if ((p10 != lb_10) && (p10 !=lb0_1) && (p10 != lb01)) lb10 = p10 ;
07134                 }
07135 
07136                         /* label point */
07137                         if (*pdest == CNCTCMP_WHITE) {
07138                                 if (lb_10) *pdest = lb_10 ;
07139                                 else if (lb0_1) *pdest = lb0_1 ;
07140                                 else if (lb01) *pdest = lb01 ;
07141                                 else if (lb10) *pdest = lb10 ;
07142                                 else {
07143                                         if (curclass < MAX_LABEL) {
07144                                                 /* avoid using rsvd label*/
07145                                                 if (curclass != CNCTCMP_WHITE) *pdest = curclass ;
07146                                                 else e_error("labeling order error at p00(%d,%d)=%d lb_10=%d, lb0_1=%d, lb01=%d, lb10=%d, cc=%d",
07147                                                                 i,j,*pdest,lb_10,lb0_1,lb01,curclass) ;
07148                                                 if (table_ok) table_ok = add_equiv_cl(cl,
07149 
07150                curclass,
07151 
07152                curclass) ;
07153                                                 curclass++ ;
07154                                         } else {
07155                                                 table_ok = FALSE ;
07156                                                 break ;
07157                                         }
07158                                 }
07159                         }
07160 
07161                         if ((*pdest != lb_10) && (*pdest !=lb0_1) &&
07162                                 (*pdest != lb01) && (*pdest !=lb10)) lb00 = *pdest ;
07163                         else lb00 = CNCTCMP_BLACK ;
07164 
07165                         /* add equivalences, avoid redundant relations*/
07166                         if (table_ok) table_ok = add_equiv_cl(cl, lb00, lb0_1) ;
07167                         if (table_ok) table_ok = add_equiv_cl(cl, lb00, lb01) ;
07168                         if (table_ok) table_ok = add_equiv_cl(cl, lb00, lb10) ;
07169                         if (table_ok) table_ok = add_equiv_cl(cl, lb00, lb_10) ;
07170                         if (table_ok) table_ok = add_equiv_cl(cl, lb_10, lb0_1) ;
07171                         if (table_ok) table_ok = add_equiv_cl(cl, lb_10, lb01) ;
07172                         if (table_ok) table_ok = add_equiv_cl(cl, lb_10, lb10) ;
07173                         if (table_ok) table_ok = add_equiv_cl(cl, lb0_1, lb01) ;
07174                         if (table_ok) table_ok = add_equiv_cl(cl, lb0_1, lb10) ;
07175                         if (table_ok) table_ok = add_equiv_cl(cl, lb01, lb10) ;
07176 
07177                         /* terminate if table saturated! */
07178                         if (!table_ok) break;
07179                 }
07180 
07181                 /* terminate if table saturated! */
07182                 if (!table_ok) break ;
07183         }
07184 
07185         /* report nb of new classes found */
07186         return curclass - first_class ;
07187 }
07188 
07189 /*--------------------------------------------------------------------------*/
07197 /*--------------------------------------------------------------------------*/
07198 static int resolve_transitivities(cleq_tab_t * cleq)
07199 {
07200         unsigned int    transitivities = 0 ;
07201         size_t                  i,
07202                                         j ;
07203 
07204         qsort(&cleq->tab[0].org,
07205                         cleq->numentries,
07206                         sizeof(cleq_t),
07207                         (qsort_fkptr_t)cleq_sort) ;
07208 
07209         for (i=0 ; i<cleq->numentries ; i++) {
07210                 for (j=i+1 ; j<cleq->numentries ; j++) {
07211                         if (cleq->tab[i].dest == cleq->tab[j].org) {
07212                                 cleq->tab[j].org = cleq->tab[i].org ;
07213                                 transitivities++ ;
07214                         }
07215                 }
07216         }
07217         return transitivities ;
07218 }
07219 /*--------------------------------------------------------------------------*/
07227 /*--------------------------------------------------------------------------*/
07228 static int eliminate_redundancies(cleq_tab_t * cleq)
07229 {
07230         unsigned int    removed = 0 ;
07231         label_pix               sorg,
07232                                         sdest ;
07233         size_t                  i,
07234                                         j ;
07235 
07236         if (!cleq->numentries) {
07237                 e_warning("eliminate_redundancies: numentries(%d)", cleq->numentries) ;
07238                 return 0 ;
07239         }
07240 
07241         for (i=0 ; i<cleq->numentries ; i++) {
07242                 sorg = cleq->tab[i].org ;
07243                 sdest= cleq->tab[i].dest ;
07244                 if (sorg && sdest) {
07245                         for (j=i+1 ; j<cleq->numentries ; j++) {
07246                                 if ((sorg == cleq->tab[j].org) && (sdest == cleq->tab[j].dest)) {
07247                                         /* clear entry */
07248                                         cleq->tab[j].org = cleq->tab[j].dest = 0 ;
07249                                         removed++ ;
07250                                 }
07251                         }
07252                 }
07253         }
07254 
07255         qsort(&cleq->tab[0].org,
07256                         cleq->numentries,
07257                         sizeof(cleq_t),
07258                         (qsort_fkptr_t)cleq_sort) ;
07259         if (removed < cleq->numentries) {
07260                 cleq->numentries -= removed ;
07261                 return removed ;
07262         } else {
07263                  e_error("eliminate_redundancies: removed (%d) larger than numentries(%d)",
07264                                         removed, cleq->numentries) ;
07265                  return -1 ;
07266         }
07267 }
07268 /*--------------------------------------------------------------------------*/
07276 /*--------------------------------------------------------------------------*/
07277 static int printcltab(cleq_tab_t *cleq)
07278 {
07279         size_t  i;
07280 
07281         e_comment(1, "Equivalences (%d entries):", cleq->numentries) ;
07282         for (i=0 ; i<cleq->numentries ; i++)
07283                 e_comment(2, "%4d\t%4d %4d", i, cleq->tab[i].org, cleq->tab[i].dest) ;
07284         return 0 ;
07285 }
07286 /*--------------------------------------------------------------------------*/
07295 /*--------------------------------------------------------------------------*/
07296 static int enforce_injectiveness(cleq_tab_t * cleq)
07297 {
07298         unsigned int    multmap = 0 ;
07299         size_t                  i,
07300                                         j ;
07301         label_pix               minlab,
07302                                         dest ;
07303 
07304         for (i=0 ; i<cleq->numentries ; i++) {
07305                 minlab = cleq->tab[i].org ;
07306                 dest = cleq->tab[i].dest ;
07307 
07308                 for (j=i+1 ; j<cleq->numentries ; j++)
07309                         if ((dest == cleq->tab[j].dest) && (cleq->tab[j].org < minlab))
07310                                 minlab = cleq->tab[j].org ;
07311 
07312                 for (j=i ; j<cleq->numentries ; j++)
07313                         if ((cleq->tab[j].org > minlab) && (dest == cleq->tab[j].dest)) {
07314                                 cleq->tab[j].dest = cleq->tab[j].org ;
07315                                 cleq->tab[j].org = minlab ;
07316                                 multmap ++ ;
07317                         }
07318         }
07319 
07320         return multmap ;
07321 }
07322 /*--------------------------------------------------------------------------*/
07331 /*--------------------------------------------------------------------------*/
07332 static long int relabel_image(
07333                 label_image_t   *       l,
07334                 cleq_tab_t              *       cl)
07335 {
07336         label_pix               maxlab = 0,
07337                                         pix ;
07338         /* correspondance table  16 -> 16 bits */
07339         ushort16        *       cw ;
07340         label_pix       *       ot ;
07341         size_t                  nextlab = FIRST_LABEL,
07342                                         max_cw_index = 0 ;
07343         size_t                  nbpix ;
07344         size_t                  i,
07345                                         j ;
07346 
07347 
07348         nbpix = l->lx * l->ly ;
07349 
07350         /* 16 bit lookup of 16 bit */
07351         cw = malloc(LUT16SIZ) ;
07352         if (cw == NULL) {
07353                 e_error("relabel_image: malloc of %d bytes failed", LUT16SIZ) ;
07354                 return -1 ;
07355         }
07356 
07357         /* mark unused entries */
07358         memset(cw, NO_CORRESPONDANCE, LUT16SIZ) ;
07359 
07360         /* Find out the index max and fill cw: cw[dest] = index */
07361         for (j=0 ; j<cl->numentries ; j++) {
07362                 if (max_cw_index < cl->tab[j].dest) max_cw_index = cl->tab[j].dest ;
07363                 cw[cl->tab[j].dest] = j ;
07364         }
07365 
07366         /* Allocate 'output table' : ot[dest] = org      */
07367         ot = calloc(CLEQSIZE, sizeof(label_pix)) ;
07368         if (ot == NULL) {
07369                 e_error("relabel_image: inable to allocate ot(%d bytes)",
07370                                   CLEQSIZE*sizeof(label_pix)) ;
07371                 free(cw) ;
07372                 return -1 ;
07373         }
07374 
07375         /* Find out maxlab and set ot : ot[i]=i for each i not in destination */
07376         /* ot[something present in destination] = 0 */
07377         for (j=0 ; j<=max_cw_index ; j++) {
07378                 if (cw[j] != NO_CORRESPONDANCE) pix = cl->tab[cw[j]].org ;
07379                 else pix = j ;
07380                ot[pix] = pix ;
07381                 if(pix > maxlab) maxlab = pix ;
07382         }
07383 
07384         /* ot[2]=2,ot[3]=0,ot[4]=4,ot[5]=5 --> ot[2]=2,ot[3]=0,ot[4]=3, ot[5]=4 */
07385         for (i=FIRST_LABEL ; i<=maxlab ; i++) if (ot[i]) ot[i] = nextlab++ ;
07386 
07387         /* Special cases */
07388         ot[(int)CNCTCMP_BLACK] = CNCTCMP_BLACK ;
07389         ot[(int)CNCTCMP_WHITE] = CNCTCMP_WHITE ;
07390 
07391         /* Update label_image object */
07392         l->nobjtot = nextlab ;
07393         l->nrealobj = l->nobjtot - FIRST_LABEL ;
07394 
07395         /* No object left in the image */
07396         if (!l->nrealobj)
07397                 e_warning("Nothing to do (found %d real objects)", l->nrealobj) ;
07398 
07399         /* Modify the input label image if necessary     */
07400         for (i=0 ; i<nbpix ; i++) {
07401                 pix = l->buf[i] ;
07402                 if (cw[pix] != NO_CORRESPONDANCE) l->buf[i] = ot[cl->tab[cw[pix]].org] ;
07403                 else if ((pix != CNCTCMP_BLACK) && (pix != CNCTCMP_WHITE))
07404                         if ((ot[pix]) && (pix != ot[pix])) {
07405                                 e_warning("relabel_image: cw[%d]=%d, ot[%d}=%d",
07406                                                 pix, cw[pix], pix, ot[pix]) ;
07407                                 pix = ot[pix] ;
07408                                 l->buf[i] = pix ;
07409                         }
07410         }
07411 
07412         /* Free and return */
07413         free(cw) ;
07414         free(ot) ;
07415         return nextlab - FIRST_LABEL ;
07416 }
07417 
07418 
07419 /*--------------------------------------------------------------------------*/
07430 /*--------------------------------------------------------------------------*/
07431 static cc_stats_t *relabel_n_calc_image(
07432                 label_image_t   *       l,
07433                 cleq_tab_t              *       cl)
07434 {
07435         label_pix               maxlab = 0,
07436                                         pix,
07437                                         maxpix ;
07438         /* correspondance table  16 -> 16 bits */
07439         ushort16        *       cw = NULL ;
07440         size_t                  nextlab = FIRST_LABEL,
07441                                         max_cw_index = 0 ;
07442         label_pix       *       ot ;
07443         cc_stats_t      *       cc = NULL,
07444                                 *       cs = NULL ;
07445         /* 4 neighbour pixels */
07446         label_pix               p0_1,
07447                                         p_10,
07448                                         p10,
07449                                         p01 ;
07450 
07451         int                             i,
07452                                         j,
07453                                         jl ;
07454 
07455         /* Test input label image */
07456         if (l == NULL) {
07457                 e_error("relabel_n_calc_image: NULL label_image provided") ;
07458                 return NULL ;
07459         }
07460 
07461         /* 16 bit lookup of 16 bit */
07462         cw = malloc(LUT16SIZ) ;
07463         if (cw == NULL) {
07464                 e_error("relabel_n_calc_image: malloc of %d bytes failed", LUT16SIZ) ;
07465                 return NULL ;
07466         }
07467 
07468         /* Mark unused entries */
07469         memset(cw, NO_CORRESPONDANCE, LUT16SIZ) ;
07470 
07471 
07472         /* Find out the index max and fill cw: cw[dest] = index */
07473         for (j=0 ; j<cl->numentries ; j++) {
07474                 if (max_cw_index < cl->tab[j].dest) max_cw_index = cl->tab[j].dest ;
07475                 cw[cl->tab[j].dest] = j ;
07476         }
07477 
07478 
07479         /* Allocate 'output table' : ot[dest] = org      */
07480         ot = calloc(CLEQSIZE, sizeof(label_pix)) ;
07481         if (ot == NULL) {
07482                 e_error("relabel_n_calc_image: unable to allocate ot(%d bytes)",
07483                                   CLEQSIZE*sizeof(label_pix)) ;
07484                 free(cw) ;
07485                 return NULL ;
07486         }
07487 
07488         /* Find out maxlab and set ot : ot[i]=i for each i not in destination */
07489         /* ot[something present in destination] = 0 */
07490         for (j=0 ; j<=max_cw_index ; j++) {
07491                 if (cw[j] != NO_CORRESPONDANCE) pix = cl->tab[cw[j]].org ;
07492                 else pix = j ;
07493                 ot[pix] = pix ;
07494                 if(pix>maxlab) maxlab = pix ;
07495         }
07496 
07497         /* ot[2]=2,ot[3]=0,ot[4]=4,ot[5]=5 --> ot[2]=2,ot[3]=0,ot[4]=3, ot[5]=4 */
07498         for(i=FIRST_LABEL ; i<=maxlab ; i++) if (ot[i]) ot[i] = nextlab++ ;
07499 
07500         /* Special cases */
07501         ot[(int)CNCTCMP_BLACK] = CNCTCMP_BLACK ;
07502         ot[(int)CNCTCMP_WHITE] = CNCTCMP_WHITE ;
07503 
07504         /* Update label_image object */
07505         l->nobjtot = nextlab ;
07506         l->nrealobj = l->nobjtot - FIRST_LABEL ;
07507 
07508         /* No object left in the image */
07509         if (!l->nrealobj)
07510                 e_warning("relabel_n_calc_image: Nothing to do (found %d real objects)",
07511                                         l->nrealobj) ;
07512 
07513         /* New statistic table */
07514         cs = new_cc_stats_table(l->nobjtot) ;
07515         if (cs == NULL) {
07516                 free(ot) ;
07517                 free(cw) ;
07518                 return NULL ;
07519         }
07520 
07521         /* Initalise features */
07522         for (i=0 ; i<l->nobjtot ; i++) {
07523                 cc = &cs[i] ;
07524                 cc ->top_y = -1 ;
07525                 cc ->rig_x = -1 ;
07526                 cc ->lef_x = l->lx ;
07527                 cc ->bottom_y = l->ly ;
07528         }
07529 
07530         /* Computing features for l->nrealobj real objects */
07531         maxpix = l->nobjtot ;
07532         for (j=0 ; j<l->ly ; j++) {
07533                 jl = j * l->lx ;
07534                 for (i=0 ; i<l->lx ; i++) {
07535                         pix = l->buf[i+jl] ;
07536                        /* Modify the input label image if necessary     */
07537                         if (cw[pix] != NO_CORRESPONDANCE) {
07538                                 l->buf[i+jl] = ot[cl->tab[cw[pix]].org] ;
07539                         } else {
07540                                 if ((pix != CNCTCMP_BLACK) && (pix != CNCTCMP_WHITE))
07541                                         if ((ot[pix]) && (pix != ot[pix])) {
07542                                                 e_warning("relabel_n_calc_image: cw[%d]=%d, ot[%d}=%d",
07543                                                                           pix, cw[pix], pix, ot[pix]) ;
07544                                                 l->buf[i+jl] = ot[pix] ;
07545                                         }
07546                         }
07547 
07548                         /* Calc features */
07549                         pix = l->buf[i+jl] ;
07550                         if (pix >= maxpix) {
07551                                  e_error("inconsistent nb of real objects=%d / label=%d",
07552                                                         l->nrealobj, pix) ;
07553                                  break ;
07554                         }
07555                         cc = &cs[pix] ;
07556 
07557                         cc->numpix++ ;
07558                         cc->gcx += i ;
07559                         cc->gcy += j ;
07560                         cc->egcx += i * i ;
07561                         cc->egcy += j * j ;
07562 
07563                         /* calculate wet surface by considering neighbour configurations*/
07564                         get4nbs(l, i, j, &p0_1, &p_10, &p01, &p10) ;
07565                         if (pix != CNCTCMP_BLACK) {
07566                                 if (p0_1 == CNCTCMP_BLACK) cc->wetsrf++ ;
07567                                 if (p_10 == CNCTCMP_BLACK) cc->wetsrf++ ;
07568                                 if (p01 == CNCTCMP_BLACK) cc->wetsrf++ ;
07569                                 if (p10 == CNCTCMP_BLACK) cc->wetsrf++ ;
07570                         } else {
07571                                 if (p0_1 != CNCTCMP_BLACK) cc->wetsrf++ ;
07572                                 if (p_10 != CNCTCMP_BLACK) cc->wetsrf++ ;
07573                                 if (p01 != CNCTCMP_BLACK) cc->wetsrf++ ;
07574                                 if (p10 != CNCTCMP_BLACK) cc->wetsrf++ ;
07575                         }
07576                         if (j < cc->bottom_y) {
07577                                 cc->bottom_y = j ;
07578                                 cc->bottom_x = i ;
07579                         }
07580                         if (j > cc->top_y) {
07581                                 cc->top_y = j ;
07582                                 cc->top_x = i ;
07583                         }
07584                         if (i > cc->rig_x) {
07585                                 cc->rig_y = j ;
07586                                 cc->rig_x = i ;
07587                         }
07588                         if (i < cc->lef_x) {
07589                                 cc->lef_y = j ;
07590                                 cc->lef_x = i ;
07591                         }
07592                }
07593         }
07594         free(cw) ;
07595         free(ot) ;
07596 
07597         /* Sum it up */
07598         for (i=0 ; i<l->nobjtot ; i++) {
07599                 cc = &cs[i] ;
07600                 if (cc->numpix) {
07601                         if (i == CNCTCMP_WHITE) {
07602                                 e_warning("object=%d has size %d pixels, image badly labeled",
07603                                                   i, cc->numpix) ;
07604                         }
07605                 } else {
07606                         if (i != CNCTCMP_WHITE) {
07607                                 e_warning("object=%d has size %d pixels",
07608                                                   i, cc->numpix) ;
07609                         }
07610                         /* don't divide by size 0*/
07611                         continue;
07612                 }
07613                 cc->gcx /= (double)cc->numpix ;
07614                 cc->gcy /= (double)cc->numpix ;
07615                 cc->egcx = sqrt(cc->egcx/(double)cc->numpix - cc->gcx) ;
07616                 cc->egcy = sqrt(cc->egcy/(double)cc->numpix - cc->gcy) ;
07617         }
07618 
07619         return cs ;
07620 }
07621 
07622 /*--------------------------------------------------------------------------*/
07634 /*--------------------------------------------------------------------------*/
07635 label_image_t * pixelmap_2_label_image(pixel_map *pm)
07636 {
07637         label_image_t   *       lab ;
07638         int                                     i ;
07639 
07640         if (pm == NULL) return NULL ;
07641 
07642         lab = new_label_image(pm->lx, pm->ly) ;
07643         for (i=0 ; i<pm->nbpix ; i++) {
07644                 lab->buf[i] = (label_pix)pm->data[i] ;
07645         }
07646         return lab ;
07647 }
07648 
07649 /*--------------------------------------------------------------------------*/
07676 /*--------------------------------------------------------------------------*/
07677 cc_stats_t * labelize_binary(label_image_t * label_image)
07678 {
07679         cc_stats_t              *       cs = NULL ;
07680         cleq_tab_t              *       cleq ;
07681         unsigned long           new_classes = 0 ;
07682         label_pix                       numobj = 0 ;
07683         int                                     Cycle = 1,
07684                                                 redundancies,
07685                                                 equival,
07686                                                 reductions ;
07687         int                                     Feature_calc_Done = FALSE ;
07688         int                                     err = 0 ;
07689 
07690    if (label_image == NULL) {
07691         e_error("NULL label_image: aborting labelize_binary") ;
07692       return NULL ;
07693    }
07694    if (label_image->buf == NULL) {
07695       e_error("NULL label_image buffer: aborting labelize_binary") ;
07696       return NULL ;
07697    }
07698 
07699         cleq = calloc(1, sizeof(cleq_tab_t)) ;
07700         if (cleq == NULL) {
07701                 e_error("labelize_binary: unable to allocate cleq (%d bytes)",
07702                                   sizeof(cleq_tab_t)) ;
07703                 return NULL ;
07704         }
07705 
07706         cleq->tab = calloc(CLEQSIZE, sizeof(cleq_t)) ;
07707         if (cleq->tab == NULL) {
07708                 e_error("labelize_binary: unable to allocate cleq->tab (%d bytes)",
07709                                   CLEQSIZE*sizeof(cleq_t)) ;
07710                 free(cleq) ;
07711                 return NULL ;
07712         }
07713 
07714         do {
07715                 new_classes = search_n_classify(label_image, FIRST_LABEL, cleq ) ;
07716                 do {
07717                         equival = resolve_transitivities(cleq) ;
07718                         redundancies = eliminate_redundancies(cleq) ;
07719                         if (redundancies < 0) {
07720                                 err = redundancies ;
07721                                 printcltab(cleq) ;
07722                                 break ;
07723                         }
07724 
07725                         reductions = enforce_injectiveness(cleq) ;
07726                 } while ((reductions || redundancies) && !err) ;
07727                 if (new_classes >= MAX_CLASSNUM) {
07728                         e_comment(2, "pass %d: merging & reordering equivalent clusters", 2) ;
07729                         numobj = relabel_image(label_image, cleq) ;
07730                         memset(cleq->tab, 0, sizeof(cleq_t)*CLEQSIZE) ;
07731                         cleq->numentries = 0 ;
07732                 } else {
07733                         cs = relabel_n_calc_image(label_image, cleq) ;
07734                         numobj = label_image->nrealobj ;
07735                         Feature_calc_Done = TRUE ;
07736                 }
07737                 if (new_classes >= MAX_CLASSNUM) e_comment(2, "cycle %d:", ++Cycle) ;
07738         } while (new_classes && new_classes >= MAX_CLASSNUM && !err && (Cycle<20)) ;
07739 
07740         if (!Feature_calc_Done) {
07741                 cs = new_cc_stats_table(numobj) ;
07742                 err = calc_ccfeatures(label_image, cs) ;
07743         }
07744         label_image->nrealobj = numobj ;
07745         free(cleq->tab) ;
07746         free(cleq) ;
07747 
07748         if (!err) return cs ;
07749         else return NULL ;
07750 }
07751 /*--------------------------------------------------------------------------*/
07761 /*--------------------------------------------------------------------------*/
07762 void free_label_image(label_image_t * to_destroy)
07763 {
07764         if (to_destroy == NULL) return ;
07765         if (to_destroy->buf != NULL) {
07766                 free(to_destroy->buf) ;
07767                 to_destroy->buf = NULL ;
07768         }
07769         free(to_destroy) ;
07770         to_destroy = NULL ;
07771         return ;
07772 }
07773 /*--------------------------------------------------------------------------*/
07784 /*--------------------------------------------------------------------------*/
07785 int dump_labelimage(
07786                 label_image_t   *       to_dump,
07787                 char                    *       filename)
07788 {
07789         FILE                    *       out ;
07790         fits_header     *       fh ;
07791         size_t                  to_write,
07792                                                 written ;
07793         char                            cval[80] ;
07794 
07795     if (filename == NULL) {
07796         e_error("no file name given: aborting dumping of label image") ;
07797         return -1 ;
07798     }
07799     if (to_dump == NULL) {
07800         e_error("NULL input label image: aborting dump") ;
07801         return -1 ;
07802     }
07803 
07804     e_comment(2, "saved %8d real objects, total labels %d",
07805                                 to_dump->nrealobj, to_dump->nobjtot) ;
07806     if (to_dump->nrealobj < 1)
07807         e_warning("Warning! No objects found in label image file %s", filename) ;
07808 
07809         fh = fits_header_default() ;
07810         sprintf(cval, "%d", BPP_16_SIGNED) ;
07811         fits_header_add(fh, "BITPIX", cval, "bits per pixel", NULL) ;
07812         sprintf(cval, "2") ;
07813         fits_header_add(fh, "NAXIS", cval, "single image", NULL) ;
07814         sprintf(cval, "%d", to_dump->lx) ;
07815         fits_header_add(fh, "NAXIS1", cval, "x axis", NULL) ;
07816         sprintf(cval, "%d", to_dump->ly) ;
07817         fits_header_add(fh, "NAXIS2", cval, "x axis", NULL) ;
07818 
07819         /* Add up number of labels */
07820         sprintf(cval, "%d", to_dump->nobjtot) ;
07821         fits_header_add(fh, "LABELS", cval, "total objects", NULL) ;
07822         /* Add number of real objects */
07823         sprintf(cval, "%d", to_dump->nrealobj) ;
07824         fits_header_add(fh, "REALOBJS", cval, "real objects", NULL) ;
07825 
07826     if ((out = fopen(filename, "a")) == NULL) {
07827         e_error("cannot append dump file [%s]", filename) ;
07828         return -1 ;
07829     }
07830         fits_dump_header(fh, out) ;
07831         fits_header_destroy(fh) ;
07832 
07833         to_write = (size_t)(to_dump->lx * to_dump->ly) ;
07834         written = fwrite(to_dump->buf,BYTESPERPIXEL(BPP_16_SIGNED), to_write, out) ;
07835         fclose(out) ;
07836 
07837     if (to_write != written) {
07838         e_error("%ld bytes written to disk instead of %ld", written, to_write) ;
07839         e_error("cannot write data to disk (disk full?): aborting") ;
07840         return -1 ;
07841     }
07842 
07843     stuff_with_zeros(filename, to_dump->lx * to_dump->ly) ;
07844     return 0 ;
07845 }
07846 /*--------------------------------------------------------------------------*/
07859 /*--------------------------------------------------------------------------*/
07860 int print_ccfeaturestable(
07861                 cc_stats_t              *       cs,
07862                 unsigned long           numentries,
07863                 FILE                    *       f,
07864                 char                    *       fin_name)
07865 {
07866         cc_stats_t              *       cc = NULL ;
07867         unsigned long           j ;
07868         signed long                     wetsrf_chksum = 0 ;
07869 
07870         if (cs == NULL) {
07871                 e_error("print_ccfeaturestable: NULL cc_stats_t provided") ;
07872                 return -1 ;
07873         }
07874         if (fin_name) fprintf(f,"# FILE %s.fits\n", fin_name) ;
07875         else fprintf(f,"# UNKNOWN input fits FILE\n") ;
07876         fprintf(f, "# Object=%d is Black (background)\n", CNCTCMP_BLACK) ;
07877         fprintf(f, "# Object=%d is White (unlabeled - all params should be zero)\n",
07878                                 CNCTCMP_WHITE);
07879         fprintf(f,
07880 "# Object  size in         Gravity           Gravity         Wet       bottom       top      left     right\n") ;
07881         fprintf(f,
07882 "# number  pixels          center            deviation       Surface\n") ;
07883         for (j=0 ; j<numentries ; j++) {
07884                 cc = &cs[j] ;
07885                 fprintf(f,
07886 "%8ld  %8d  %8.2f,%8.2f  %8.2f,%8.2f  %8ld %4d,%4d %4d,%4d %4d,%4d %4d,%4d\n",
07887                                 j,
07888                                 cc->numpix,
07889                                 (double)cc->gcx,
07890                                 (double)cc->gcy,
07891                                 (double)cc->egcx,
07892                                 (double)cc->egcy,
07893                                 cc->wetsrf,
07894                                 (int)cc->bottom_x,
07895                                 (int)cc->bottom_y,
07896                                 (int)cc->top_x,
07897                                 (int)cc->top_y,
07898                                 (int)cc->lef_x,
07899                                 (int)cc->lef_y,
07900                                 (int)cc->rig_x,
07901                                 (int)cc->rig_y) ;
07902                 if (j != CNCTCMP_BLACK) wetsrf_chksum += cc->wetsrf ;
07903         }
07904         fprintf(f, "# Total %ld Objects including Black and White\n", numentries) ;
07905         fprintf(f, "# WetSurface BLACK - total is %ld\n",
07906                           wetsrf_chksum - cs[(int)CNCTCMP_BLACK].wetsrf) ;
07907         return 0 ;
07908 }
07909 
07910 
07911 
07912 
07913 /*---------------------------------------------------------------------------*/
07932 /*---------------------------------------------------------------------------*/
07933 static int select_valid_arcs(
07934                 int                                     min_arc_length,
07935                 int                                     max_arc_width,
07936                 cc_stats_t              *       cc_stats,
07937                 int                             nobjs,
07938                 int                     *       validarcs,
07939                 int                     **      lut_cc_2_arc,
07940                 int                     **      lut_arc_2_cc,
07941                 orientation_t           arc_orientation)
07942 {
07943         int             *       tl = NULL ;
07944         cc_stats_t      *       c ;
07945         int                     ymin,
07946                                         ymax ;
07947         double                  arclen,
07948                                         archeight ;
07949         int                             err = 0 ;
07950         int                             i,
07951                                         j ;
07952 
07953         /* Initialize    */
07954         ymin = 0 ;
07955         ymax = 999999999 ;
07956         *lut_cc_2_arc = NULL ;
07957         *lut_arc_2_cc = NULL ;
07958 
07959         /**************/
07960         /* DEBUG CASE */
07961         if (debug_active() >= ARC_DEBUG_LEVEL)
07962                 e_comment(0, "select_valid_arcs:  min_arc_length=%d max_arc_width=%d",
07963                                         min_arc_length, max_arc_width) ;
07964         /**************/
07965 
07966         /**************/
07967         /* DEBUG CASE */
07968         if ((ymax-ymin<min_arc_length) && (debug_active()>=ARC_DEBUG_LEVEL))
07969                 e_warning("valid arcs: arc range (%d) < min_range (%d)",
07970                                 ymax-ymin,
07971                                 min_arc_length) ;
07972         /**************/
07973 
07974         /**************/
07975         /* DEBUG CASE */
07976         if (debug_active()>=ARC_DEBUG_LEVEL)
07977                 e_comment(1,"select_valid_arcs: %d valid arcs found, limlo=%d limhi=%d",
07978                         *validarcs, ymin, ymax) ;
07979         /**************/
07980 
07981 
07982         /* Allocate look-up table between objects and arcs */
07983         *lut_cc_2_arc = malloc(nobjs*sizeof(int)) ;
07984         j = 0 ;
07985         *validarcs = 0 ;
07986 
07987         /* For each object, test if it is big enough to be a valid arc */
07988         /* Count the number of valid arcs and update lut_cc_2_arc array */
07989         for (i=FIRST_LABEL ; i<nobjs ; i++) {
07990                 c = &(cc_stats[i]) ;
07991                 switch(arc_orientation) {
07992                         case VERTICAL:
07993                         arclen = c->top_y - c->bottom_y + 1 ;
07994                         archeight = c->rig_x - c->lef_x + 1 ;
07995                         if ((arclen > min_arc_length) && (archeight < max_arc_width) &&
07996 
07997                (c->lef_x >0)) {
07998                                 (*lut_cc_2_arc)[i] = j++ ;
07999                                 (*validarcs)++ ;
08000                                 if (ymin < c->bottom_y) ymin = c->bottom_y ;
08001                                 if (ymax > c->top_y) ymax = c->top_y ;
08002                         } else {
08003                                 (*lut_cc_2_arc)[i] = -1 ;
08004                         }
08005                         break ;
08006 
08007                         case HORIZONTAL:
08008                         arclen = c->rig_x - c->lef_x + 1 ;
08009                         archeight = c->top_y - c->bottom_y + 1 ;
08010                         if ((arclen > min_arc_length) && (archeight < max_arc_width) &&
08011 
08012                (c->bottom_y>0)) {
08013                                 (*lut_cc_2_arc)[i] = j++ ;
08014                                 (*validarcs)++ ;
08015                                 if (ymin < c->lef_x) ymin = c->lef_x ;
08016                                 if (ymax > c->rig_x) ymax = c->rig_x ;
08017                         } else {
08018                                 (*lut_cc_2_arc)[i] = -1 ;
08019                         }
08020                         break ;
08021 
08022                         default:
08023                         e_error("select_valid_arcs: unsupported orientation %d",
08024                                         arc_orientation) ;
08025                         return -1 ;
08026                 }
08027         }
08028 
08029         /* No valid arc found */
08030         if (!(*validarcs)) return err ;
08031 
08032         /* Allocate look-up table between arcs and objects */
08033         *lut_arc_2_cc = calloc(*validarcs, sizeof(int)) ;
08034 
08035         /* Associate to each arc, its object number */
08036         tl = *lut_arc_2_cc ;
08037         for (i=FIRST_LABEL ; i<nobjs ; i++) {
08038                 c = &(cc_stats[i]) ;
08039                 switch(arc_orientation) {
08040                         case VERTICAL:
08041                         arclen = c->top_y - c->bottom_y + 1 ;
08042                         archeight = c->rig_x - c->lef_x + 1 ;
08043                         if ((arclen>min_arc_length) && (archeight<max_arc_width) &&
08044                                                                                                 (c->lef_x>0)) *tl++ = i ;
08045                         break ;
08046 
08047                         case HORIZONTAL:
08048                         arclen = c->rig_x - c->lef_x + 1 ;
08049                         archeight = c->top_y - c->bottom_y + 1 ;
08050                         if ((arclen>min_arc_length) && (archeight<max_arc_width) &&
08051                                                                                                 (c->bottom_y >0)) *tl++ = i ;
08052                         break ;
08053                 }
08054         }
08055 
08056         return err ;
08057 }
08058 /*--------------------------------------------------------------------------*/
08068 /*--------------------------------------------------------------------------*/
08069 int select_objs(
08070                 label_image_t   *       l,
08071                 int                     *       objlist,
08072                 int                             nobjs)
08073 {
08074         int                     err = 0 ;
08075         long int                i,
08076                                         nbpix ;
08077         int                     *       lut = NULL ;
08078 
08079         if (l == NULL) {
08080                 e_error("select_objs : label image NULL") ;
08081                 return -1 ;
08082         }
08083         if (objlist == NULL) {
08084                 e_error("select_objs : objlist NULL") ;
08085                 return -1 ;
08086         }
08087         if (nobjs > l->nrealobj) {
08088                 e_error("select_objs : objects know %d objects to select %d",
08089                                 l->nrealobj, nobjs) ;
08090                 return -1 ;
08091         }
08092         nbpix = l->lx * l->ly ;
08093         if (!nbpix) {
08094                 e_error("select_objs : nbpix %d", nbpix) ;
08095                 return -1 ;
08096         }
08097         lut = calloc(l->nobjtot, sizeof(int)) ;
08098         if (lut == NULL) {
08099                 e_error("select_objs : failed allocating lut for %d objects", nobjs) ;
08100                 return -1 ;
08101         }
08102         for (i=0 ; i<nobjs ; i++) lut[objlist[i]] = objlist[i] ;
08103 
08104         e_comment(0, "select_objs : deleting %d objects", l->nrealobj - nobjs) ;
08105         l->nrealobj = nobjs ;
08106         for (i=0;i<nbpix;i++) l->buf[i] = lut[l->buf[i]] ;
08107         free(lut) ;
08108 
08109         return err ;
08110 }
08111 
08112 /*---------------------------------------------------------------------------*/
08137 /*---------------------------------------------------------------------------*/
08138 static cc_stats_t       *       detect_arcs(
08139         arc_param_t             *       arc_param,
08140         OneImage                *       in,
08141         int                             *       n_arcs,
08142         int                     **      lut_arc_2_cc,
08143         int                     **      lut_cc_2_arc,
08144         label_image_t   **      label_image)
08145 {
08146         OneImage        *       filt_img  = NULL,
08147                                 *       collapsed = NULL ;
08148         pixel_map       *       thresh    = NULL ;
08149         FILE            *       txtfile   = NULL ;
08150         cc_stats_t      *       cc_stats  = NULL ;
08151         pixelvalue      threshold,
08152                                         fillval ;
08153         double                  median_val, med_val2,
08154                                         sigma ;
08155         int                             collapse_direction ;
08156         int                             err = 0 ;
08157 
08158         /* Default values for output parameters */
08159         *lut_arc_2_cc = NULL ;
08160         *lut_cc_2_arc = NULL ;
08161         *label_image  = NULL ;
08162 
08163 
08164         /* Set up the collapse direction */
08165         switch(arc_param->arc_orientation){
08166                 case VERTICAL:
08167                 collapse_direction = 0 ;
08168                 break ;
08169 
08170                case HORIZONTAL:
08171                 collapse_direction = 1 ;
08172                 break ;
08173                 default:
08174                 e_error("detect_arcs: unknown orientation %d",
08175                                 arc_param->arc_orientation) ;
08176                 return NULL ;
08177         }
08178 
08179         /* Clear zones to be ignored (to avoid false detections) */
08180         median_val = image_median_stat(in, &sigma) ;
08181         fillval = median_val-sigma/2.0 ;
08182 
08183         /**************/
08184         /* DEBUG CASE */
08185         if (debug_active() >= ARC_DEBUG_LEVEL)
08186                 e_error("detect_arcs:  fillval=%f median=%f sigma=%f",
08187                                 (float)fillval, (float)median_val, (float)sigma) ;
08188         /**************/
08189 
08190         if (arc_param->badtop > 0)
08191                 fillrect_in_image(in, fillval,
08192                                                         0, in->lx-1, 0, arc_param->badtop-1) ;
08193         if (arc_param->badbot < in->ly -1)
08194                 fillrect_in_image(in, fillval,
08195                                                         0, in->lx-1, arc_param->badbot+1, in->ly-1) ;
08196         if (arc_param->badleft > 0 )
08197                 fillrect_in_image(in, fillval,
08198                                                         0, arc_param->badleft-1, 0, in->ly-1) ;
08199         if (arc_param->badright< in->lx -1)
08200                 fillrect_in_image(in, fillval,
08201                                                         arc_param->badright+1,in->lx-1, 0, in->ly-1) ;
08202 
08203 
08204         /* Median filter in detection direction */
08205         switch(arc_param->arc_orientation) {
08206                 case VERTICAL:
08207                 filt_img = image_filter_vertical_median(in, arc_param->median_size) ;
08208                 break ;
08209                 case HORIZONTAL:
08210                 filt_img = image_filter_horizontal_median(in, arc_param->median_size) ;
08211                 break ;
08212         }
08213 
08214         /* Subtract a low-pass */
08215         /*if (subtract_lowpass(arc_param->arc_orientation,
08216                                                         arc_param->window_size,
08217                                                         filt_img) == -1) {
08218                 destroy_image(filt_img) ;
08219                 return NULL ;
08220         }*/
08221 
08222         /* Get relevant stats for thresholding */
08223     median_val = image_median_stat(filt_img, &sigma) ;
08224        /**************/
08225         /* DEBUG CASE */
08226         if (debug_active() >= ARC_DEBUG_LEVEL)
08227                 e_error("detect_arcs:  median=%f, sigma=%f",
08228                                 (float)median_val, (float)sigma) ;
08229         /**************/
08230 
08231         /* Collapse the image */
08232         collapsed = collapse_image_median(filt_img, collapse_direction, 0, 0) ;
08233         med_val2 = image_median_stat(collapsed, &sigma) ;
08234 
08235         /* Correct median_val and sigma if necessary */
08236         if (median_val < TRESH_MEDIAN_MIN) median_val = TRESH_MEDIAN_MIN ;
08237         if (sigma > TRESH_SIGMA_MAX) sigma = TRESH_SIGMA_MAX ;
08238 
08239         /* Set the threshold */
08240         threshold = median_val + sigma*arc_param->thresh_fact;
08241 
08242         /**************/
08243         /* DEBUG CASE */
08244         if (debug_active() >= ARC_DEBUG_LEVEL)
08245                 e_error("using threshold: thresh=%f med_val2=%f", threshold, med_val2*2) ;
08246         /**************/
08247 
08248         /* Threshold to keep only the arcs - use of the collapsed image */
08249         err = threshold_one_dim(filt_img,
08250                                                                 med_val2*2,
08251                                                                 collapsed->data,
08252                                                                 arc_param->arc_orientation,
08253                                                                 0.0) ;
08254         destroy_image(collapsed);
08255         if (err) {
08256                 destroy_image(filt_img);
08257                 return NULL;
08258         }
08259 
08260         /* Binarize the image */
08261         thresh = thresh_image_to_pixelmap(filt_img, threshold, MAX_PIX_VALUE) ;
08262         if (thresh==NULL) {
08263                 destroy_image(filt_img) ;
08264                 return NULL ;
08265         }
08266 
08267         /* Test if there are enough good pixels */
08268         if (thresh->ngoodpix < arc_param->min_good_pixels) {
08269                 e_error("detect_arcs: too few (%d) white pixels", thresh->ngoodpix) ;
08270 
08271                 /**************/
08272                 /* DEBUG CASE */
08273                 if (debug_active() >= ARC_DEBUG_LEVEL)
08274                         save_image_to_fits(filt_img, "med.fits", ARC_PARAM_PTYPE) ;
08275                 /**************/
08276 
08277                 destroy_image(filt_img) ;
08278                 destroy_pixelmap(thresh) ;
08279                 return NULL ;
08280         }
08281         destroy_image(filt_img) ;
08282 
08283         /* Convert the pixel map to a label image */
08284         *label_image = pixelmap_2_label_image(thresh) ;
08285         destroy_pixelmap(thresh) ;
08286 
08287         /* Identify the different clusters and give each a label */
08288         cc_stats = labelize_binary(*label_image) ;
08289         if (cc_stats==NULL) {
08290                 if (*label_image) free_label_image(*label_image) ;
08291                 return NULL ;
08292         }
08293 
08294         /**************/
08295         /* DEBUG CASE */
08296         if (debug_active() >= ARC_DEBUG_LEVEL) {
08297                 dump_labelimage(*label_image, "lab.fits") ;
08298                 txtfile = fopen("lab.txt", "w") ;
08299                 print_ccfeaturestable(cc_stats,
08300                                                                 (*label_image)->nobjtot,
08301                                 txtfile,
08302                                                                 "detect_arcs.dummy") ;
08303                 fclose(txtfile) ;
08304         }
08305         /**************/
08306 
08307         /* Select only relevant arcs, create corresponding LUT's */
08308         printf("arc_param->min_arclen=%f arc_param->max_arcwidth=%f, (*label_image)->nobjtot=%d,n_arcs=%d, arc_param->arc_orientation=%d\n",
08309                 arc_param->min_arclen,arc_param->max_arcwidth, (*label_image)->nobjtot, n_arcs, arc_param->arc_orientation);
08310         err = select_valid_arcs(arc_param->min_arclen,
08311                                                         arc_param->max_arcwidth,
08312                                                         cc_stats,
08313                                                         (*label_image)->nobjtot,
08314                                                         n_arcs,
08315                                                         lut_cc_2_arc,
08316                                                         lut_arc_2_cc,
08317                                                         arc_param->arc_orientation) ;
08318         if (err) {
08319                 if (*lut_cc_2_arc!=NULL) free(*lut_cc_2_arc) ;
08320                 *lut_cc_2_arc = NULL ;
08321         if (*lut_arc_2_cc!=NULL) free(*lut_arc_2_cc) ;
08322                 *lut_arc_2_cc = NULL ;
08323                 free(cc_stats) ;
08324                 return NULL ;
08325         }
08326 
08327         if (*n_arcs < arc_param->min_nb_arcs) {
08328                 e_error("detect_arcs: number of detected arcs %d is too low <%d",
08329                                 *n_arcs, arc_param->min_nb_arcs) ;
08330 
08331                 /**************/
08332                 /* DEBUG CASE */
08333                 if (debug_active() >= ARC_DEBUG_LEVEL && *n_arcs > 0) {
08334                         select_objs(*label_image, *lut_arc_2_cc, *n_arcs) ;
08335                         dump_labelimage(*label_image, "sel.fits") ;
08336                 }
08337                 /**************/
08338 
08339                 if (*lut_cc_2_arc != NULL) free(*lut_cc_2_arc) ;
08340                 *lut_cc_2_arc = NULL ;
08341         if (*lut_arc_2_cc != NULL) free(*lut_arc_2_cc) ;
08342                 *lut_arc_2_cc = NULL ;
08343                 free(cc_stats) ;
08344                 return NULL ;
08345         }
08346         return cc_stats ;
08347 }
08348 
08349 /*---------------------------------------------------------------------------*/
08359 /*---------------------------------------------------------------------------*/
08360 double * dist_engine( arc_param_t *arc_param, OneImage *org, float offset)
08361 {
08362         dpoint                  **      grid2d       = NULL ;
08363         double3                 *       lamgrid      = NULL ;
08364         cc_stats_t              *       cc_stats     = NULL ;
08365         label_image_t   *       label_image  = NULL ;
08366         OneImage                *       in           = NULL,
08367                                         *       masked       = NULL,
08368                                         *       collapsed    = NULL ;
08369         double                  *       refgrid      = NULL ;
08370         double                  *       refgrid_tmp  = NULL ;
08371         double                          refgrid_total;
08372         double                  *       distpoly     = NULL ;
08373         double                          mean_squared_error ;
08374         double          pixval ;
08375         int                                     n_arcs ;
08376         int                             *       lut_arc_2_cc = NULL,
08377                                         *       lut_cc_2_arc = NULL ;
08378         int                                     ncoeffs ;
08379         int                                     min_arc_range ;
08380         int                                     left,
08381                                                 right,
08382                                                 top,
08383                                                 bottom ;
08384         int                 err = 0 ;
08385         int                 n_lines = 0 ;
08386 
08387         int                             i,
08388                                                 j ;
08389 
08390         /* Test if input image is valid  */
08391         if (org==NULL){
08392                 e_error("dist_engine: input image NULL");
08393                 return NULL;
08394         }
08395         in = copy_image(org);
08396 
08397         /* Detect the arcs in the input frame */
08398         cc_stats = detect_arcs(arc_param,
08399                                                         in,
08400                                                         &n_arcs,
08401                                                         &lut_arc_2_cc,
08402                                                         &lut_cc_2_arc,
08403                                                         &label_image) ;
08404 
08405 
08406 
08407 
08408         if (n_arcs<32){
08409                 e_error("wrong nummber of arcs: %d",n_arcs);
08410                 return NULL;
08411         }
08412         if (cc_stats==NULL) {
08413                 if (label_image!=NULL) free_label_image(label_image);
08414                 destroy_image(in);
08415                 return NULL;
08416         }
08417         if ((lut_arc_2_cc==NULL) || (lut_cc_2_arc==NULL)) {
08418                 free(cc_stats);
08419                 free_label_image(label_image);
08420                 destroy_image(in);
08421                 return NULL;
08422         }
08423         if (n_arcs<=0) {
08424                 free(cc_stats);
08425                 free_label_image(label_image);
08426                 destroy_image(in);
08427                 free(lut_cc_2_arc);
08428                 return NULL;
08429         }
08430 
08431         /**************/
08432         /* DEBUG CASE */
08433         if (debug_active()>=ARC_DEBUG_LEVEL) {
08434                 select_objs(label_image, lut_arc_2_cc, n_arcs);
08435                 dump_labelimage(label_image, "sel.fits");
08436         }
08437         /**************/
08438         /* Find out the limit objects positions */
08439         err = get_extreme_obj_coor(cc_stats,
08440                                                                 label_image->nobjtot,
08441                                                                 lut_arc_2_cc,
08442                                                                 n_arcs,
08443                                                                 &left,
08444                                                                 &right,
08445                                                                 &top,
08446                                                                 &bottom) ;
08447 
08448         /* Abort if the detected arcs are too concentrated in the same zone */
08449         switch(arc_param->arc_orientation) {
08450 
08451                 case VERTICAL:
08452                 min_arc_range = in->lx/ARC_RANGE_FACT ;
08453                 if (right-left<min_arc_range) {
08454                         e_error("dist_engine: too narrow arc range (%d-%d)<%d",
08455                                         right, left, min_arc_range) ;
08456                         free(cc_stats) ;
08457                         free_label_image(label_image) ;
08458                         destroy_image(in) ;
08459                         free(lut_cc_2_arc) ;
08460                         free(lut_arc_2_cc) ;
08461                         return NULL ;
08462                 }
08463                 break ;
08464                case HORIZONTAL:
08465                 min_arc_range = in->ly/ARC_RANGE_FACT ;
08466                 if (top-bottom<min_arc_range) {
08467                         e_error("dist_engine: too narrow arc range (%d-%d)<%d",
08468                                         bottom, top, min_arc_range) ;
08469                         free(cc_stats) ;
08470                         free_label_image(label_image) ;
08471                         destroy_image(in) ;
08472                         free(lut_cc_2_arc) ;
08473                         free(lut_arc_2_cc) ;
08474                         return NULL ;
08475                 }
08476         }
08477 
08478         /* Create a 2-D deformation grid with detected arcs */
08479         grid2d = get_positions(arc_param,
08480                                                         in,
08481                                                         n_arcs,
08482                                                         lut_cc_2_arc,
08483                                                         label_image,
08484                                                         cc_stats) ;
08485         free(lut_cc_2_arc);
08486         lut_cc_2_arc=NULL;
08487         if (grid2d==NULL){
08488                 free_label_image(label_image);
08489                 destroy_image(in);
08490                 return NULL;
08491         }
08492 
08493         /* lamgrid is (x,y,X) where (x,y) is the arc coord., X is the ideal pos */
08494         lamgrid = double3_new(arc_param->n_calib*n_arcs);
08495         /* refgrid will contain the arc positions, ie the Xs of lamgrid */
08496         refgrid = calloc(n_arcs, sizeof(double));
08497         refgrid_tmp = calloc(n_arcs, sizeof(double));
08498 
08499         /**************/
08500         /* DEBUG CASE */
08501         if (debug_active()>=ARC_DEBUG_LEVEL)
08502                 e_comment(0,"creating surface grid for %d arcs",n_arcs);
08503         /**************/
08504 
08505         switch(arc_param->grid_ref_point) {
08506 
08507                 case ARC_GRID_REF_1DMAX:
08508                 e_comment(0,"using collapsed peaks for arc reference grid");
08509                 masked = copy_image(in);
08510 
08511                 switch(arc_param->arc_orientation) {
08512 
08513                         case VERTICAL:
08514                         for (i=0 ; i<n_arcs ; i++) {
08515                                 /* Keep only the current arc - mask the others */
08516                                 err = mask_objs(in,
08517                                                                 masked,
08518                                                                 label_image,
08519                                                                 &(lut_arc_2_cc[i]),
08520                                                                 1,
08521 
08522                                                                 0.0) ;
08523 
08524                                 /**************/
08525                                 /* DEBUG CASE */
08526                                 if (debug_active() >= ARC_DEBUG_LEVEL && !i)
08527                                         save_image_to_fits(masked,"masked.fits",ARC_PARAM_PTYPE);
08528                                 /**************/
08529 
08530                                 collapsed = collapse_image(masked, 0);
08531 
08532                                 /* First estimation of the current arc position */
08533                                 refgrid[i] = cc_stats[lut_arc_2_cc[i]].gcx ;
08534 
08535                                 /* Refining of the current arc position */
08536                                 if ((refgrid[i] > 5) && (refgrid[i] < in->lx-5)) {
08537                                         refgrid[i] += function1d_find_centroid(
08538                                                                                 &(collapsed->data[(int)refgrid[i]-5]),
08539                                                                                 11) - 5 ;
08540                                         /* e_comment(0, "gcx = %f max = %f", */
08541                                                                         /* (float) cc_stats[lut_arc_2_cc[i]].gcx, */
08542                                                                         /* (float)refgrid[i]) ;  */
08543                                 } else {
08544                                         e_warning("arc %d (%f) too close to border for function1d_find_centroid()",
08545                                                                 i,
08546                                                                 (float)refgrid[i]) ;
08547                                 }
08548                                 destroy_image(collapsed) ;
08549                         }
08550                         destroy_image(masked) ;
08551                         break;
08552 
08553                         case HORIZONTAL:
08554                         for (i=0 ; i<n_arcs ; i++) {
08555                                 /* Keep only the current arc - mask the others */
08556                                 err = mask_objs(in,
08557                                                                 masked,
08558                                                                 label_image,
08559                                                                 &(lut_arc_2_cc[i]),
08560                                                                 1,
08561                                                                 0.0) ;
08562 
08563                                 /**************/
08564                                 /* DEBUG CASE */
08565                                 if (debug_active() >= ARC_DEBUG_LEVEL && !i)
08566                                         save_image_to_fits(masked,"masked.fits",ARC_PARAM_PTYPE);
08567                                 /**************/
08568 
08569                                 collapsed = collapse_image(masked, 1) ;
08570 
08571                                 /* First estimation of the current arc position */
08572                                 refgrid[i] = cc_stats[lut_arc_2_cc[i]].gcy ;
08573 
08574                                 /* Refining of the current arc position */
08575                                 if ((refgrid[i] > 5) && (refgrid[i] < in->ly-5)){
08576                                        refgrid[i] += function1d_find_centroid(
08577                                                                                 &(collapsed->data[(int)refgrid[i]-5]),
08578                                                                                 11) - 5 ;
08579                                         /* e_comment(0,"gcy = %f max = %f", */
08580                                                                         /* (float) cc_stats[lut_arc_2_cc[i]].gcy, */
08581                                                                         /* (float)refgrid[i]);  */
08582                                 } else {
08583                                         e_warning("arc %d (%f)too close to border for function1d_find_centroid()",
08584                                                                 i,
08585                                                                 (float)refgrid[i]) ;
08586                                 }
08587                                 destroy_image(collapsed) ;
08588                         }
08589                         destroy_image(masked) ;
08590                         break ;
08591                 }
08592                 break;
08593 
08594                 case ARC_GRID_REF_GRAV_CENT:
08595                 e_comment(0, "using gravity centers for arc reference grid") ;
08596                 for (i=0 ; i<n_arcs ; i++) refgrid[i]=cc_stats[lut_arc_2_cc[i]].gcx ;
08597                 break;
08598 
08599                 case ARC_FIND_MAX_CENT:
08600                 masked = copy_image(in);
08601                 refgrid_total=0;
08602                 for(i=0 ;i<n_arcs;i++) {
08603                         /* Keep only the current arc - mask the others */
08604                         err = mask_objs(in,
08605                                                         masked,
08606                                                         label_image,
08607                                                         &(lut_arc_2_cc[i]),
08608                                                         1,
08609                                                         0.0) ;
08610 
08611                         /**************/
08612                         /* DEBUG CASE */
08613                         if (debug_active() >= ARC_DEBUG_LEVEL && !i)
08614                                 save_image_to_fits(masked,"masked.fits",ARC_PARAM_PTYPE);
08615                         /**************/
08616 
08617                         /*collapsed = collapse_image(masked, 0);*/
08618 
08619 
08620                         collapsed = new_image(masked->lx, 1) ;
08621                         for (j=0 ; j<collapsed->lx ; j++) {
08622                                 pixval=get_sum_pixels_on_vig(masked,
08623                                             j+1,
08624                                             (masked->lx/2)-50,
08625                                             j+1,
08626                                             (masked->lx/2)+50) ;
08627                                 collapsed->data[j] = (pixelvalue)pixval ;
08628                         }
08629                         /* First estimation of the current arc position */
08630                         if ((refgrid[i] > 5) && (refgrid[i] < in->lx-5)) {
08631                                 refgrid[i] += function1d_find_centroid(
08632                                                                         &(collapsed->data[(int)refgrid[i]-5]),
08633                                                                         11) - 5 ;
08634                                 /* e_comment(0, "gcx = %f max = %f", */
08635                                                                 /* (float) cc_stats[lut_arc_2_cc[i]].gcx, */
08636                                                                 /* (float)refgrid[i]) ;  */
08637                         } else {
08638                                 e_warning("arc %d (%f) too close to border for function1d_find_centroid()",
08639                                                         i,
08640                                                         (float)refgrid[i]) ;
08641                         }
08642                         /*refgrid_total += refgrid[i]-refgrid[0];*/
08643                         destroy_image(collapsed) ;
08644                 }
08645                 for(i=0 ;i<n_arcs;i++)
08646                    refgrid_tmp[i]=refgrid[i];
08647                 qsort(refgrid_tmp, n_arcs, sizeof(double), comp_q);
08648                 /*for(i=4 ;i<n_arcs;i=i+3)
08649                 {
08650                    refgrid_total += refgrid_tmp[i]-refgrid_tmp[1];
08651                 }
08652                 for(i=0 ;i<n_arcs;i++) {
08653                   refgrid[i]=refgrid[i]*((masked->lx/32.0)*(31*32/2)/refgrid_total);
08654                   e_error("%g",refgrid[i]);
08655                 }*/
08656                 n_lines = n_arcs/32.0;
08657                 for(i=n_lines ;i<n_arcs;i=i+n_lines)
08658                 {
08659                    for(j=0;j<n_lines;j++){
08660                      refgrid_total += refgrid_tmp[i+j]-refgrid_tmp[j];
08661                    }
08662                 }
08663                 for(i=0 ;i<n_arcs;i++) {
08664                   refgrid[i]=refgrid[i]*((masked->lx/32.0)*n_lines*(31*32/2))/refgrid_total-offset;
08665                   refgrid_tmp[i]=refgrid_tmp[i]*((masked->lx/32.0)*n_lines*(31*32/2))/refgrid_total-offset;
08666                   /*if(i>0)
08667                     e_error("%g",refgrid_tmp[i]-refgrid_tmp[i-1]);*/
08668                 }
08669                 destroy_image(masked) ;
08670                 break;
08671 
08672                 default:
08673                 e_error("unknown grid reference point type %d",
08674                                 arc_param->grid_ref_point) ;
08675                 free(refgrid);
08676                 free(refgrid_tmp);
08677                 double3_del(lamgrid);
08678                 destroy_image(in);
08679                 free_label_image(label_image);
08680                 return NULL;
08681         }
08682         free_label_image(label_image);
08683 
08684 
08685         /* Fill lamgrid with refgrid and grid2d */
08686         switch(arc_param->arc_orientation) {
08687 
08688                 case VERTICAL:
08689                 for (i=0 ; i<n_arcs ; i++) {
08690 
08691                         /**************/
08692                         /* DEBUG CASE */
08693                         if (debug_active()>=ARC_DEBUG_LEVEL)
08694                                 printf("\nArc %3d (%4.1f): ", i, (float)refgrid[i]) ;
08695                         /**************/
08696 
08697                         for (j=0 ; j<arc_param->n_calib ; j++) {
08698                                 lamgrid->x[i*arc_param->n_calib+j] = refgrid[i] ;
08699                                 lamgrid->y[i*arc_param->n_calib+j] = grid2d[i][j].y ;
08700                                 lamgrid->z[i*arc_param->n_calib+j] = grid2d[i][j].x ;
08701 
08702                                 /**************/
08703                                 /* DEBUG CASE */
08704                                 if (debug_active()>=ARC_DEBUG_LEVEL)
08705                                         printf("%4.0f %6.2f   ",
08706                                                 lamgrid->y[i*arc_param->n_calib+j],
08707                                                 lamgrid->z[i*arc_param->n_calib+j]) ;
08708                                 /**************/
08709 
08710                         }
08711                         free(grid2d[i]) ;
08712                 }
08713                 break ;
08714 
08715 
08716                 case HORIZONTAL:
08717                 for (i=0 ; i<n_arcs ; i++) {
08718 
08719                         /**************/
08720                         /* DEBUG CASE */
08721                         if (debug_active()>=ARC_DEBUG_LEVEL)
08722                                 printf("\nArc %d (%4.1f): ", i, (float)refgrid[i]) ;
08723                         /**************/
08724 
08725                         for (j=0 ; j<arc_param->n_calib ; j++) {
08726                                 lamgrid->y[i*arc_param->n_calib+j] = refgrid[i] ;
08727                                 lamgrid->x[i*arc_param->n_calib+j] = grid2d[i][j].x ;
08728                                 lamgrid->z[i*arc_param->n_calib+j] = grid2d[i][j].y ;
08729 
08730                                 /**************/
08731                                 /* DEBUG CASE */
08732                                 if (debug_active()>=ARC_DEBUG_LEVEL)
08733                                         printf("%6.2f %4.0f  ",
08734                                                 lamgrid->x[i*arc_param->n_calib+j],
08735                                                 lamgrid->z[i*arc_param->n_calib+j]) ;
08736                                 /**************/
08737 
08738                         }
08739                         free(grid2d[i]) ;
08740                 }
08741                 break ;
08742         }
08743         free(grid2d);
08744         free(refgrid);
08745         free(refgrid_tmp);
08746 
08747         /* Compute the 2D polynomial fit */
08748         distpoly = fit_surface_polynomial(lamgrid,
08749                                                                                 "(0,0) (1,0) (0,1) (1,1) (2,0) (0,2) (2,1) (1,2)",
08750                                                                                 3,
08751                                                                                 &ncoeffs,
08752                                                                                 &mean_squared_error) ;
08753 
08754         /****************/
08755         /* VERBOSE CASE */
08756         /*if (verbose_active())*/
08757                 arc_print_results(distpoly, mean_squared_error, n_arcs) ;
08758         /****************/
08759 
08760         /* Free and return the 2D polynomial */
08761         double3_del(lamgrid) ;
08762         free(cc_stats) ;
08763         free(lut_arc_2_cc) ;
08764         lut_arc_2_cc = NULL ;
08765         destroy_image(in) ;
08766         return distpoly ;
08767 }
08768 
08769 /*---------------------------------------------------------------------------*/
08788 /*---------------------------------------------------------------------------*/
08789 poly2d * compute_distortion(
08790         OneImage                *       in,
08791         int                                     badleft,
08792         int                                     badright,
08793         int                                     badtop,
08794         int                                     badbot,
08795         int             orientation,
08796         float           offset)
08797 {
08798         double          *       distpoly  = NULL ;
08799         arc_param_t     *       arc_param = NULL ;
08800         poly2d          *       poly_u    = NULL ;
08801         orientation_t           arc_orientation;
08802 
08803         if(orientation==0)
08804                 arc_orientation = HORIZONTAL;
08805         else
08806                 arc_orientation = VERTICAL;
08807 
08808         /* Set default values in the arc_param object */
08809         /*arc_param= setup_default_arc_param(   ARC_PARAM_WINDOWSIZE,
08810                                                                                 ARC_PARAM_MEDIAN_SIZE,
08811                                                                                 ARC_PARAM_THRESHFACT,
08812                                                                                 ARC_GRID_REF_GRAV_CENT) ;*/
08813         arc_param= setup_default_arc_param(     ARC_PARAM_WINDOWSIZE,
08814                                                                                 ARC_PARAM_MEDIAN_SIZE,
08815                                                                                 ARC_PARAM_THRESHFACT,
08816                                                                                 ARC_FIND_MAX_CENT) ;
08817         if (arc_param==NULL) return NULL;
08818 
08819         arc_param->badtop   = badtop ;
08820         arc_param->badbot   = badbot ;
08821         arc_param->badleft  = badleft ;
08822         arc_param->badright = badright ;
08823         arc_param->arc_orientation = arc_orientation ;
08824         arc_param->max_arcwidth = ARC_PARAM_MAXARCWIDTH ;
08825 
08826         switch(arc_param->arc_orientation) {
08827 
08828                 case HORIZONTAL:
08829                 arc_param->min_arclen = in->lx/ARC_PARAM_MINARCLENFACT ;
08830                 break ;
08831 
08832                 case VERTICAL:
08833                 arc_param->min_arclen = in->ly/ARC_PARAM_MINARCLENFACT ;
08834                 break ;
08835 
08836                 default:
08837                 e_error("unknown arc orientation %d,(%d=vertical,%d=horizontal",
08838                                 arc_param->arc_orientation, VERTICAL, HORIZONTAL) ;
08839                 return NULL ;
08840         }
08841 
08842         /* Apply the distortion computation engine */
08843         distpoly = dist_engine(arc_param, in, offset) ;
08844         free(arc_param) ;
08845         if (distpoly==NULL) return NULL ;
08846 
08847         /* Construct fill and return the output object (poly2d) */
08848         poly_u=malloc(sizeof(poly2d)) ;
08849         poly_u->nc = 8 ;
08850         poly_u->px=calloc(8, sizeof(int)) ;
08851         poly_u->py=calloc(8, sizeof(int)) ;
08852         poly_u->px[0] = 0 ;
08853         poly_u->px[1] = 1 ;
08854         poly_u->px[2] = 0 ;
08855         poly_u->px[3] = 1 ;
08856         poly_u->px[4] = 2 ;
08857         poly_u->px[5] = 0 ;
08858         poly_u->px[6] = 2 ;
08859         poly_u->px[7] = 1 ;
08860         poly_u->py[0] = 0 ;
08861         poly_u->py[1] = 0 ;
08862         poly_u->py[2] = 1 ;
08863         poly_u->py[3] = 1 ;
08864         poly_u->py[4] = 0 ;
08865         poly_u->py[5] = 2 ;
08866         poly_u->py[6] = 1 ;
08867         poly_u->py[7] = 2 ;
08868         poly_u->c = distpoly ;
08869 
08870         return poly_u ;
08871 }
08872 
08873 
08874 
08875 /*-------------------------------------------------------------------------*/
08886 /*--------------------------------------------------------------------------*/
08887 
08888 
08889 int is_ascii_list(char * name)
08890 {
08891         FILE    *       in ;
08892         char            line[FILENAMESZ] ;
08893         int                     i ;
08894 
08895         /* Get the first FILENAMESZ characters from the file */
08896         in = fopen(name, "r") ;
08897         if (in == NULL)
08898                 return -1 ;
08899         fread(line, sizeof(char), FILENAMESZ-1, in) ;
08900         fclose(in) ;
08901         /* Check out that they are all ASCII characters */
08902         for (i=0 ; i<FILENAMESZ ; i++) {
08903                 if ((i>1) && (line[i]=='\n')) return 1 ;
08904                 if (!isascii(line[i])) return 0 ;
08905         }
08906         return 1 ;
08907 }
08908 
08909 /*-------------------------------------------------------------------------*/
08923 /*--------------------------------------------------------------------------*/
08924 
08925 OneCube *
08926 load_cube_ASCII_list(char * listname)
08927 {
08928         OneCube         *       loaded ;
08929         int                             nfiles ;
08930         char            **      filenames ;
08931 
08932 
08933         /*
08934          * Load the list of all file names contained in the list
08935          */
08936         filenames = get_ASCII_list_names(listname, &nfiles);
08937         if ((nfiles<1) || (filenames==NULL)) {
08938                 e_error("reading ASCII list [%s]: aborting load", listname);
08939                 return NULL ;
08940         }
08941 
08942         /*
08943          * Outsource the loading to load_cube_string_list()
08944          */
08945         loaded = load_cube_string_list(filenames, nfiles) ;
08946         free_ASCII_list_names(filenames, nfiles);
08947         return loaded ;
08948 }
08949 
08950 
08951 
08952 /*---------------------------------------------------------------------------*/
08967 /*---------------------------------------------------------------------------*/
08968 arc_param_t     *       setup_default_arc_param(
08969                 int                     window_size,
08970                 int             median_size,
08971                 double          thresh_fact,
08972                 grid_refpoint_t grid_ref_point)
08973 {
08974         arc_param_t     *       arc_param = NULL ;
08975 
08976         /* Allocate object */
08977         arc_param = calloc(1, sizeof(arc_param_t)) ;
08978         if (arc_param==NULL) {
08979                 e_error("setup_default_arc_param: arc_param NULL") ;
08980                 return NULL ;
08981         }
08982 
08983         /* Set values */
08984         arc_param->median_size          = median_size ;
08985         arc_param->window_size          = window_size ;
08986         arc_param->grid_ref_point       = grid_ref_point ;
08987 
08988         /* Set defaults for other values */
08989         arc_param->n_calib                      = ARC_PARAM_NBSAMPLES ;
08990         arc_param->thresh_fact          = thresh_fact ;
08991         arc_param->min_nb_arcs          = ARC_PARAM_MINNBARCS ;
08992         arc_param->min_good_pixels      = ARC_PARAM_MINGOODPIX ;
08993 
08994         /* Return */
08995         return arc_param ;
08996 }
08997 

Generated on Wed Oct 26 13:08:53 2005 for SINFONI Pipeline Reference Manual by doxygen1.2.13.1 written by Dimitri van Heesch, © 1997-2001