00001
00002
00003
00004
00005
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
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
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
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
00090
00091
00092
00093
00094 static size_t emem_ramlimit = 1024 ;
00095 static size_t emem_vm_limit = 1024 ;
00096 static size_t emem_vm_pagesize = 1024 ;
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
00105
00106
00107
00108
00109
00110
00111
00112
00113
00114
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
00129 static malloc_cell * malloc_table = NULL ;
00130
00131
00132
00133
00134
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
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
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
00265 if (file_exists(filename) != 1) {
00266 e_error("file %s not found", filename) ;
00267 return((cube_info*)NULL) ;
00268 }
00269
00270
00271 fileinfo = malloc(sizeof(cube_info)) ;
00272
00273
00274 fileinfo->headersize = fits_get_header_size(filename);
00275 if (fileinfo->headersize<1) {
00276 free(fileinfo);
00277 return NULL ;
00278 }
00279
00280
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
00294
00295
00296
00297
00298
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
00306
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
00325
00326
00327 if (naxes == 1) {
00328 fileinfo->ly = 1 ;
00329 } else {
00330
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
00354
00355 if (naxes <= 2) {
00356 fileinfo->n_im = 1 ;
00357 } else {
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
00374 }
00375
00376
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
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
00400 if (atoi(value)==0) {
00401 e_warning("BSCALE=0 in header: using 1");
00402 fileinfo->b_scale = 1.0 ;
00403 }
00404 }
00405
00406
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
00414
00415
00416
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
00491
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
00561 test1 = where[len] ;
00562 test2 = where[len+1] ;
00563
00564 if (test1=='=') break ;
00565
00566
00567
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
00576 mmap_close(mm);
00577 return NULL ;
00578 }
00579
00580 where += 80 ;
00581 bufcount +=80;
00582 if (bufcount>=fileinfo.st_size){
00583
00584
00585 mmap_close(mm);
00586 return NULL ;
00587 }
00588
00589 }
00590
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
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
00688
00689
00690 if (!strncmp(line, "END ", 4)) {
00691 return NULL ;
00692 }
00693
00694
00695
00696
00697 memset(value, 0, 81);
00698
00699 if (!strncmp(line, "HISTORY ", 8) || !strncmp(line, " ", 8)) {
00700 i=7 ;
00701
00702 while (line[i]==' ' && i<80) i++ ;
00703 if (i>=80) return NULL ;
00704 from=i ;
00705
00706
00707 to=79 ;
00708 while (line[to]==' ') to-- ;
00709
00710 strncpy(value, line+from, to-from+1);
00711
00712 value[to-from+1] = (char)0;
00713 return value ;
00714 } else if (!strncmp(line, "COMMENT ", 8)) {
00715
00716
00717
00718
00719 i=7 ;
00720 while (line[i]==' ' && i<80) i++ ;
00721 if (i>=80) return NULL ;
00722 from=i ;
00723
00724
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
00735 strncpy(value, line+from, to-from+1);
00736
00737 value[to-from+1] = (char)0;
00738 return value ;
00739 }
00740
00741
00742
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
00764
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
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
00795 strncpy(value, line+from, to-from+1);
00796
00797 value[to-from+1] = (char)0;
00798
00799
00800
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
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
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
00882
00883
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
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
00929
00930
00931
00932
00933
00934
00935
00936
00937
00938
00939
00940
00941
00942
00943
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957
00958
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
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
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
01086 #if (BYTE_ORDERING_SCHEME == INTEL)
01087
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
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
01106 #if (BYTE_ORDERING_SCHEME == INTEL)
01107
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
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
01127 #if (BYTE_ORDERING_SCHEME == INTEL)
01128
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
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
01151
01152
01153
01154
01155 if (sizeof(pixelvalue)<sizeof(double)) {
01156 e_warning("loss of precision in input") ;
01157 }
01158 #if (BYTE_ORDERING_SCHEME == INTEL)
01159
01160
01161
01162
01163
01164
01165 for (i=0 ; i<nbpix ; i++) {
01166
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
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
01183 for (i=0 ; i<nbpix ; i++) {
01184
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
01237 if (list==NULL || np<1) return NULL ;
01238
01239
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
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
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
01417 if (d->data != NULL)
01418 free(d->data) ;
01419 } else {
01420
01421
01422
01423
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
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
01472 if (d->plane != NULL)
01473 free(d->plane) ;
01474
01475
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
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
01522 if (incube==NULL) return NULL ;
01523 if (incube->np<3) {
01524
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
01805 if (s==NULL) return NULL ;
01806
01807 pretty[0] = (char)0 ;
01808 if (s[0]!='\'') return s ;
01809
01810
01811 i=1 ;
01812 j=0 ;
01813
01814 while (s[i]==' ') {
01815 if (i==(int)strlen(s)) break ;
01816 i++ ;
01817 }
01818 if (i>=(int)(strlen(s)-1)) return pretty ;
01819
01820 while (i<(int)strlen(s)) {
01821 if (s[i]=='\'') {
01822 i++ ;
01823 }
01824 pretty[j]=s[i];
01825 i++ ;
01826 j++ ;
01827 }
01828
01829 pretty[j+1]=(char)0;
01830
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
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
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
01929
01930
01931 for (plane=0 ; plane<in_cube->np ; plane++) {
01932 timeline[plane] = in_cube->plane[plane]->data[pos] ;
01933 }
01934
01935 pixel_qsort(timeline, in_cube->np) ;
01936 acc_val = 0.0 ;
01937
01938
01939
01940
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
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
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
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
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
02231 image_out = new_image(image1->lx, image1->ly) ;
02232
02233
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
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
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
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
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
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
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
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
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
02817
02818 mX = least_sq_mx(mA,mB) ;
02819
02820
02821
02822 close_mx(mA) ;
02823 close_mx(mB) ;
02824
02825
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
02844
02845 if (mse != NULL) {
02846 err = 0.00 ;
02847 for (i=0 ; i<np ; i++) {
02848 y = c[0] ;
02849
02850
02851
02852 for (k=1 ; k<=poly_deg ; k++) {
02853 xp = ipow(list[i].x, k) ;
02854 y += c[k] * xp ;
02855 }
02856
02857
02858
02859 xp = ipow(list[i].y - y, 2) ;
02860 err += xp ;
02861 }
02862
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
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
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
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
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
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
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
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
03378
03379
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
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
03416
03417
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
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
03476 if (!strncmp(line, " ", 8)) {
03477 strcpy(key, " ");
03478 return key ;
03479 }
03480
03481
03482
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
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
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
03519 strncpy(key, line, i) ;
03520
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
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
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
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
03593 while (line[i]==' ') i++ ;
03594 from=i ;
03595
03596
03597
03598
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
03610 strncpy(comment, line+from, to-from+1);
03611
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
03676 if (!is_blank_line(line)) {
03677
03678
03679 key = fits_getkey(line);
03680 val = fits_getvalue(line);
03681 com = fits_getcomment(line);
03682
03683
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
03693 fits_header_append(hdr, key, val, com, line);
03694
03695 if (strlen(key)==3)
03696 if (key[0]=='E' &&
03697 key[1]=='N' &&
03698 key[2]=='D')
03699 break ;
03700 }
03701
03702 where += 80 ;
03703
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
03770
03771
03772
03773
03774
03775
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
03789 k = malloc(sizeof(keytuple));
03790
03791 k->key = strdup2(key) ;
03792
03793 k->val = NULL ;
03794 if (val!=NULL) {
03795 if (strlen(val)>0) {
03796 k->val = strdup2(val);
03797 }
03798 }
03799
03800 k->com = NULL ;
03801 if (com!=NULL) {
03802 if (strlen(com)>0) {
03803 k->com = strdup2(com) ;
03804 }
03805 }
03806
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
03817
03818
03819
03820
03821
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
03837
03838
03839
03840
03841
03842
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
03857
03858
03859
03860
03861
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
03937 if (d->plane != NULL)
03938 free(d->plane) ;
03939
03940
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
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
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
04018
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
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
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
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
04205 for (j=0 ; j < ly_out ; j++) {
04206 for (i=0 ; i< lx_out ; i++) {
04207
04208
04209 x = poly2d_compute(poly_u, (double)i, (double)j);
04210 y = poly2d_compute(poly_v, (double)i, (double)j);
04211
04212
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
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
04229 tabx = (x - (double)px) * (double)(TABSPERPIX) ;
04230 taby = (y - (double)py) * (double)(TABSPERPIX) ;
04231
04232
04233
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
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
04266 image_out->data[i+j*lx_out] = (pixelvalue)(cur/sumrs) ;
04267
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 ;
04514 inv_np = 1.00 / (double)np ;
04515
04516
04517
04518
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
04533
04534 reverse_tanh_kernel(x, np) ;
04535
04536
04537
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
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
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
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
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
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
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
04901
04902
04903
04904
04905
04906
04907
04908
04909
04910
04911 if ((ptype == BPP_IEEE_FLOAT) &&
04912 (sizeof(pixelvalue)*BITSPERBYTE == 32)) {
04913
04914
04915
04916
04917
04918
04919
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
04940
04941
04942
04943
04944
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
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
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
05119 if (image_in==NULL) return NULL ;
05120
05121
05122 if ((fabs(shift_x)<1e-2) && (fabs(shift_y)<1e-2))
05123 return copy_image(image_in) ;
05124
05125
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
05155
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
05163
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
05177
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
05195
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
05203
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
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
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
05913 if (appended==NULL || filename==NULL) return -1 ;
05914
05915
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
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
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
05977 if (size_before==0) {
05978 e_error("NULL file size: aborting FITS file zero padding") ;
05979 return ;
05980 }
05981
05982
05983 if ((size_before % (size_t)FITS_BLOCK_SIZE) == 0)
05984 return ;
05985
05986
05987 added_zeros = (size_t)FITS_BLOCK_SIZE -
05988 (size_before % (size_t)FITS_BLOCK_SIZE) ;
05989
05990 zeros = malloc(added_zeros * sizeof(uchar8)) ;
05991 memset(zeros, 0, added_zeros) ;
05992
05993
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
06005 written = fwrite(zeros, sizeof(uchar8), added_zeros, stuffed) ;
06006
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, ®ex_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, ®ex_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, ®ex_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
06158 memset(line, ' ', 80);
06159 if (key==NULL) return ;
06160
06161
06162 if (!strcmp(key, "END")) {
06163
06164 sprintf(line, "END") ;
06165 return ;
06166 }
06167
06168 if (!strcmp(key, "HISTORY") ||
06169 !strcmp(key, "COMMENT") ||
06170 !strncmp(key, " ", 8)) {
06171
06172 sprintf(line, "%s ", key);
06173 if (val==NULL)
06174 return ;
06175
06176
06177 len = strlen(val);
06178
06179
06180
06181
06182 if (strlen(val)>72) {
06183 len=72 ;
06184 }
06185 strncpy(line+8, val, len);
06186 return ;
06187 }
06188
06189
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
06198 if (com==NULL) {
06199 strcpy(ccom, "no comment");
06200 } else {
06201 strcpy(ccom, com);
06202 }
06203
06204
06205 if (!strncmp(key, "HIERARCH", 8)) {
06206 hierarch ++ ;
06207 }
06208
06209
06210
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
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
06238
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
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
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
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
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
06360 for (i=0 ; i<nbpix ; i++) {
06361 if (p_source[i] > (pixelvalue)UCHAR_MAX) {
06362
06363 p_dest[i] = (BYTE)UCHAR_MAX ;
06364 } else {
06365 if (p_source[i] < (pixelvalue)0) {
06366
06367 mem_dest[i] = (BYTE)0 ;
06368 } else {
06369
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
06380 shortintpix = (short16)(SHRT_MAX);
06381 } else {
06382 if (p_source[i] < (pixelvalue)SHRT_MIN) {
06383
06384 shortintpix = (short16)SHRT_MIN ;
06385 } else {
06386
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
06403 longintpix = (long32)LONG_MAX ;
06404 } else {
06405 if (p_source[i] < (pixelvalue)LONG_MIN) {
06406
06407 longintpix = (long32)LONG_MIN ;
06408 } else {
06409
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
06427
06428
06429
06430
06431
06432
06433
06434
06435 floatpix = (float)p_source[i] ;
06436 longintpix = *((long*)&floatpix) ;
06437 #if (BYTE_ORDERING_SCHEME == INTEL)
06438
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
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
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
06521 fh = fits_header_default();
06522
06523 sprintf(cval, "%d", to_save->orig_ptype);
06524 fits_header_add(fh, "BITPIX", cval, "Bits per pixel", NULL);
06525
06526
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
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
06543
06544
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
06549 fits_header_add(fh, "ECLIPSE", "1", "created by eclipse", NULL);
06550 fits_header_add(fh, "ORIGIN", "eclipse", "created by eclipse", NULL);
06551
06552 write_history_in_header(fh, to_save->history, to_save->n_comments) ;
06553
06554
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
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
06580
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
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
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
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
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
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
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
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
07003 nbpix = lx * ly ;
07004
07005
07006 lab = (label_image_t*)calloc(1, sizeof(label_image_t)) ;
07007
07008
07009 lab->lx = lx ;
07010 lab->ly = ly ;
07011 lab->nrealobj = 0 ;
07012 lab->nobjtot = 0 ;
07013
07014
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
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
07048 size_t i,
07049 j ;
07050
07051 size_t jl ;
07052
07053 int table_ok ;
07054
07055 label_pix curclass ;
07056
07057 label_pix * pdest ;
07058
07059 label_pix p0_1,
07060 p_10,
07061 p10,
07062 p01 ;
07063
07064 label_pix lb0_1,
07065 lb_10,
07066 lb00,
07067 lb10,
07068 lb01 ;
07069
07070
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
07090 switch(p_10) {
07091 case CNCTCMP_BLACK:
07092
07093 break ;
07094
07095 case CNCTCMP_WHITE:
07096
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
07107 break ;
07108
07109 case CNCTCMP_WHITE:
07110
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
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
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
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
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
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
07178 if (!table_ok) break;
07179 }
07180
07181
07182 if (!table_ok) break ;
07183 }
07184
07185
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
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
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
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
07358 memset(cw, NO_CORRESPONDANCE, LUT16SIZ) ;
07359
07360
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
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
07376
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
07385 for (i=FIRST_LABEL ; i<=maxlab ; i++) if (ot[i]) ot[i] = nextlab++ ;
07386
07387
07388 ot[(int)CNCTCMP_BLACK] = CNCTCMP_BLACK ;
07389 ot[(int)CNCTCMP_WHITE] = CNCTCMP_WHITE ;
07390
07391
07392 l->nobjtot = nextlab ;
07393 l->nrealobj = l->nobjtot - FIRST_LABEL ;
07394
07395
07396 if (!l->nrealobj)
07397 e_warning("Nothing to do (found %d real objects)", l->nrealobj) ;
07398
07399
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
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
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
07446 label_pix p0_1,
07447 p_10,
07448 p10,
07449 p01 ;
07450
07451 int i,
07452 j,
07453 jl ;
07454
07455
07456 if (l == NULL) {
07457 e_error("relabel_n_calc_image: NULL label_image provided") ;
07458 return NULL ;
07459 }
07460
07461
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
07469 memset(cw, NO_CORRESPONDANCE, LUT16SIZ) ;
07470
07471
07472
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
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
07489
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
07498 for(i=FIRST_LABEL ; i<=maxlab ; i++) if (ot[i]) ot[i] = nextlab++ ;
07499
07500
07501 ot[(int)CNCTCMP_BLACK] = CNCTCMP_BLACK ;
07502 ot[(int)CNCTCMP_WHITE] = CNCTCMP_WHITE ;
07503
07504
07505 l->nobjtot = nextlab ;
07506 l->nrealobj = l->nobjtot - FIRST_LABEL ;
07507
07508
07509 if (!l->nrealobj)
07510 e_warning("relabel_n_calc_image: Nothing to do (found %d real objects)",
07511 l->nrealobj) ;
07512
07513
07514 cs = new_cc_stats_table(l->nobjtot) ;
07515 if (cs == NULL) {
07516 free(ot) ;
07517 free(cw) ;
07518 return NULL ;
07519 }
07520
07521
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
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
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
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
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
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
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
07820 sprintf(cval, "%d", to_dump->nobjtot) ;
07821 fits_header_add(fh, "LABELS", cval, "total objects", NULL) ;
07822
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
07954 ymin = 0 ;
07955 ymax = 999999999 ;
07956 *lut_cc_2_arc = NULL ;
07957 *lut_arc_2_cc = NULL ;
07958
07959
07960
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
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
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
07983 *lut_cc_2_arc = malloc(nobjs*sizeof(int)) ;
07984 j = 0 ;
07985 *validarcs = 0 ;
07986
07987
07988
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
08030 if (!(*validarcs)) return err ;
08031
08032
08033 *lut_arc_2_cc = calloc(*validarcs, sizeof(int)) ;
08034
08035
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
08159 *lut_arc_2_cc = NULL ;
08160 *lut_cc_2_arc = NULL ;
08161 *label_image = NULL ;
08162
08163
08164
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
08180 median_val = image_median_stat(in, &sigma) ;
08181 fillval = median_val-sigma/2.0 ;
08182
08183
08184
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
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
08215
08216
08217
08218
08219
08220
08221
08222
08223 median_val = image_median_stat(filt_img, &sigma) ;
08224
08225
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
08232 collapsed = collapse_image_median(filt_img, collapse_direction, 0, 0) ;
08233 med_val2 = image_median_stat(collapsed, &sigma) ;
08234
08235
08236 if (median_val < TRESH_MEDIAN_MIN) median_val = TRESH_MEDIAN_MIN ;
08237 if (sigma > TRESH_SIGMA_MAX) sigma = TRESH_SIGMA_MAX ;
08238
08239
08240 threshold = median_val + sigma*arc_param->thresh_fact;
08241
08242
08243
08244 if (debug_active() >= ARC_DEBUG_LEVEL)
08245 e_error("using threshold: thresh=%f med_val2=%f", threshold, med_val2*2) ;
08246
08247
08248
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
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
08268 if (thresh->ngoodpix < arc_param->min_good_pixels) {
08269 e_error("detect_arcs: too few (%d) white pixels", thresh->ngoodpix) ;
08270
08271
08272
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
08284 *label_image = pixelmap_2_label_image(thresh) ;
08285 destroy_pixelmap(thresh) ;
08286
08287
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
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
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
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
08391 if (org==NULL){
08392 e_error("dist_engine: input image NULL");
08393 return NULL;
08394 }
08395 in = copy_image(org);
08396
08397
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
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
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
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
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
08494 lamgrid = double3_new(arc_param->n_calib*n_arcs);
08495
08496 refgrid = calloc(n_arcs, sizeof(double));
08497 refgrid_tmp = calloc(n_arcs, sizeof(double));
08498
08499
08500
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
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
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
08533 refgrid[i] = cc_stats[lut_arc_2_cc[i]].gcx ;
08534
08535
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
08541
08542
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
08556 err = mask_objs(in,
08557 masked,
08558 label_image,
08559 &(lut_arc_2_cc[i]),
08560 1,
08561 0.0) ;
08562
08563
08564
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
08572 refgrid[i] = cc_stats[lut_arc_2_cc[i]].gcy ;
08573
08574
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
08580
08581
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
08604 err = mask_objs(in,
08605 masked,
08606 label_image,
08607 &(lut_arc_2_cc[i]),
08608 1,
08609 0.0) ;
08610
08611
08612
08613 if (debug_active() >= ARC_DEBUG_LEVEL && !i)
08614 save_image_to_fits(masked,"masked.fits",ARC_PARAM_PTYPE);
08615
08616
08617
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
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
08635
08636
08637 } else {
08638 e_warning("arc %d (%f) too close to border for function1d_find_centroid()",
08639 i,
08640 (float)refgrid[i]) ;
08641 }
08642
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
08649
08650
08651
08652
08653
08654
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
08667
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
08686 switch(arc_param->arc_orientation) {
08687
08688 case VERTICAL:
08689 for (i=0 ; i<n_arcs ; i++) {
08690
08691
08692
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
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
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
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
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
08756
08757 arc_print_results(distpoly, mean_squared_error, n_arcs) ;
08758
08759
08760
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
08809
08810
08811
08812
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
08843 distpoly = dist_engine(arc_param, in, offset) ;
08844 free(arc_param) ;
08845 if (distpoly==NULL) return NULL ;
08846
08847
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
08896 in = fopen(name, "r") ;
08897 if (in == NULL)
08898 return -1 ;
08899 fread(line, sizeof(char), FILENAMESZ-1, in) ;
08900 fclose(in) ;
08901
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
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
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
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
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
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
08995 return arc_param ;
08996 }
08997