00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019 #include "eclipse.h"
00020 #include "cnct_comp.h"
00021
00022
00023
00024
00025
00026 #define CLEQSIZE 65534
00027 #define MAXUI16 65535
00028
00029
00030 #define LUT16SIZ ((MAXUI16+1)*sizeof(ushort16))
00031 #define NO_CORRESPONDANCE (0xFFFF)
00032
00033
00034
00035
00036
00037 typedef struct _CLASS_EQUIVALENCE_ {
00038 label_pix org ;
00039 label_pix dest;
00040 } cleq_t;
00041
00042
00043 typedef struct _CLASS_EQUIVALENCE_TABLE_ {
00044 size_t numentries ;
00045 cleq_t *tab;
00046 } cleq_tab_t;
00047
00048 typedef int(*qsort_fkptr_t)(const void*, const void*) ;
00049
00050
00051
00052
00053
00054 static inline void get4nbs( label_image_t *, size_t, size_t, label_pix *,
00055 label_pix *, label_pix *, label_pix *) ;
00056
00057 static inline void get8nbs(label_image_t *, size_t, size_t, label_pix *,
00058 label_pix *, label_pix *, label_pix *, label_pix *, label_pix *,
00059 label_pix *, label_pix *) ;
00060
00061 static cc_stats_t *new_cc_stats_table(size_t) ;
00062 static int add_equiv_cl(cleq_tab_t *, label_pix, label_pix) ;
00063 static int cleq_sort(const cleq_t *, const cleq_t *) ;
00064 static int resolve_transitivities(cleq_tab_t *) ;
00065 static int eliminate_redundancies(cleq_tab_t *) ;
00066 static int enforce_injectiveness(cleq_tab_t *) ;
00067 static int printcltab(cleq_tab_t *) ;
00068 static long int relabel_image(label_image_t *, cleq_tab_t *) ;
00069 static cc_stats_t * relabel_n_calc_image(label_image_t *, cleq_tab_t *) ;
00070 static int calc_ccfeatures( label_image_t *, cc_stats_t *) ;
00071 static label_pix search_n_classify(label_image_t *, label_pix, cleq_tab_t *) ;
00072 static int get_fits_info(char *, size_t *, size_t *, int *, size_t *) ;
00073
00074
00075
00076
00077
00078
00079
00095
00096 static inline void get4nbs(
00097 label_image_t * l,
00098 size_t i,
00099 size_t j,
00100 label_pix * p0_1,
00101 label_pix * p_10,
00102 label_pix * p01,
00103 label_pix * p10)
00104 {
00105 size_t center ;
00106
00107 center = i + j * l->lx ;
00108
00109 if (j > 0) *p0_1 = l->buf[center - l->lx] ;
00110 else *p0_1 = CNCTCMP_BLACK ;
00111
00112 if (i > 0) *p_10 = l->buf[center - 1] ;
00113 else *p_10 = CNCTCMP_BLACK ;
00114
00115 if (i < l->lx-1) *p10 = l->buf[center + 1] ;
00116 else *p10 = CNCTCMP_BLACK ;
00117
00118 if (j < l->ly-1) *p01 = l->buf[center + l->lx] ;
00119 else *p01 = CNCTCMP_BLACK ;
00120 }
00121
00122
00141
00142 static inline void get8nbs(
00143 label_image_t * l,
00144 size_t i,
00145 size_t j,
00146 label_pix * p_1_1,
00147 label_pix * p0_1,
00148 label_pix * p1_1,
00149 label_pix * p_10,
00150 label_pix * p01,
00151 label_pix * p_11,
00152 label_pix * p10,
00153 label_pix * p11)
00154 {
00155 size_t center ;
00156 size_t cl ;
00157
00158 center = i + j * l->lx ;
00159
00160 if (j > 0) {
00161 cl = center - l->lx ;
00162 if (i > 0) *p_1_1 = l->buf[cl - 1] ;
00163 else *p_1_1 = CNCTCMP_BLACK ;
00164 *p0_1 = l->buf[cl] ;
00165 if (i < l->lx-1) *p1_1 = l->buf[cl + 1] ;
00166 else *p1_1 = CNCTCMP_BLACK ;
00167 } else *p_1_1 = *p0_1 = *p1_1 = CNCTCMP_BLACK ;
00168
00169 if (i > 0) *p_10 = l->buf[center-1] ;
00170 else *p_10 = CNCTCMP_BLACK ;
00171
00172 if (i < l->lx-1) *p10 = l->buf[center + 1] ;
00173 else *p10 = CNCTCMP_BLACK ;
00174
00175 if (j < l->ly-1) {
00176 cl = center + l->lx ;
00177 if (i>0) *p_11 = l->buf[cl - 1] ;
00178 else *p_11 = CNCTCMP_BLACK ;
00179 *p01 = l->buf[cl] ;
00180 if (i < l->lx-1) *p11 = l->buf[cl + 1] ;
00181 else *p11 = CNCTCMP_BLACK ;
00182 } else *p_11 = *p01 = *p11 = CNCTCMP_BLACK ;
00183 }
00184
00185
00186
00195
00196 label_image_t * new_label_image(
00197 size_t lx,
00198 size_t ly)
00199 {
00200 label_image_t * lab ;
00201 size_t nbpix ;
00202
00203
00204 if ((lx <= 0) || (lx > MAX_COLUMN_NUMBER) ||
00205 (ly <= 0) || (ly > MAX_LINE_NUMBER)) {
00206 e_error("size out of bounds: cannot create label image") ;
00207 e_error("requested size is [%d x %d]", lx, ly) ;
00208 return NULL ;
00209 }
00210
00211
00212 nbpix = lx * ly ;
00213
00214
00215 lab = (label_image_t*)calloc(1, sizeof(label_image_t)) ;
00216
00217
00218 lab->lx = lx ;
00219 lab->ly = ly ;
00220 lab->nrealobj = 0 ;
00221 lab->nobjtot = 0 ;
00222
00223
00224 lab->buf = malloc(nbpix * sizeof(label_pix)) ;
00225 if (lab == NULL) {
00226 e_error("memory allocation failure: aborting label image creation") ;
00227 free(lab) ;
00228 return NULL ;
00229 }
00230
00231
00232 return lab ;
00233 }
00234
00235
00236
00246
00247 static cc_stats_t * new_cc_stats_table(size_t numobj)
00248 {
00249 cc_stats_t * cc_stats_table = NULL ;
00250
00251
00252 if (!numobj) e_warning("new_cc_stats: nothing to do (numobj=%d)", numobj) ;
00253
00254
00255 cc_stats_table = calloc(numobj+2, sizeof(cc_stats_t)) ;
00256 if (cc_stats_table == NULL) e_error("new_cc_stats: mem full (numobj=%d)",
00257 numobj) ;
00258
00259
00260 return cc_stats_table ;
00261 }
00262
00263
00264
00265
00275
00276 void free_label_image(label_image_t * to_destroy)
00277 {
00278 if (to_destroy == NULL) return ;
00279 if (to_destroy->buf != NULL) {
00280 free(to_destroy->buf) ;
00281 to_destroy->buf = NULL ;
00282 }
00283 free(to_destroy) ;
00284 to_destroy = NULL ;
00285 return ;
00286 }
00287
00288
00289
00299
00300 static int add_equiv_cl(
00301 cleq_tab_t * cl,
00302 label_pix lab0,
00303 label_pix lab1)
00304 {
00305 label_pix sorg,
00306 sdest ;
00307 size_t i = 0 ;
00308
00309 if (!lab0) return TRUE ;
00310 if (!lab1) return TRUE ;
00311 if ((lab0 == CNCTCMP_WHITE) || (lab1 == CNCTCMP_WHITE)) {
00312 e_error("Internal add_equiv_cl error : attempt to use label WHITE(%d)",
00313 CNCTCMP_WHITE) ;
00314 return FALSE ;
00315 }
00316
00317 if (lab0 < lab1) {
00318 sorg = lab0 ;
00319 sdest = lab1 ;
00320 } else {
00321 sorg = lab1 ;
00322 sdest = lab0 ;
00323 }
00324
00325
00326 while ((i < CLEQSIZE) && cl->tab[i].org) {
00327
00328 if ((cl->tab[i].org == sorg) && (cl->tab[i].dest == sdest)) return TRUE ;
00329 else i++ ;
00330 }
00331
00332
00333 if (i < CLEQSIZE) {
00334
00335 if (i != cl->numentries) {
00336 e_error("internal add_equiv_cl error") ;
00337 e_error("inconsistent number of entries(%d), expected(%d)",
00338 cl->numentries, i) ;
00339 } else {
00340 cl->tab[i].org = sorg ;
00341 cl->tab[i].dest = sdest ;
00342 cl->numentries++ ;
00343 return TRUE ;
00344 }
00345 }
00346
00347
00348 return FALSE ;
00349 }
00350
00351
00352
00353
00354
00355
00356 static int cleq_sort(
00357 const cleq_t * c1,
00358 const cleq_t * c2)
00359 {
00360
00361 if (c1->org && c2->org) {
00362 if (c1->org > c2->org) return 1 ;
00363 else if (c1->org < c2->org) return -1 ;
00364 else return 0 ;
00365 }
00366
00367
00368 if (c2->org) return 1 ;
00369
00370 if (c1->org) return -1 ;
00371
00372 return 0 ;
00373 }
00374
00375
00376
00384
00385 static int resolve_transitivities(cleq_tab_t * cleq)
00386 {
00387 unsigned int transitivities = 0 ;
00388 size_t i,
00389 j ;
00390
00391 qsort(&cleq->tab[0].org,
00392 cleq->numentries,
00393 sizeof(cleq_t),
00394 (qsort_fkptr_t)cleq_sort) ;
00395
00396 for (i=0 ; i<cleq->numentries ; i++) {
00397 for (j=i+1 ; j<cleq->numentries ; j++) {
00398 if (cleq->tab[i].dest == cleq->tab[j].org) {
00399 cleq->tab[j].org = cleq->tab[i].org ;
00400 transitivities++ ;
00401 }
00402 }
00403 }
00404 return transitivities ;
00405 }
00406
00407
00408
00416
00417 static int eliminate_redundancies(cleq_tab_t * cleq)
00418 {
00419 unsigned int removed = 0 ;
00420 label_pix sorg,
00421 sdest ;
00422 size_t i,
00423 j ;
00424
00425 if (!cleq->numentries) {
00426 e_warning("eliminate_redundancies: numentries(%d)", cleq->numentries) ;
00427 return 0 ;
00428 }
00429
00430 for (i=0 ; i<cleq->numentries ; i++) {
00431 sorg = cleq->tab[i].org ;
00432 sdest= cleq->tab[i].dest ;
00433 if (sorg && sdest) {
00434 for (j=i+1 ; j<cleq->numentries ; j++) {
00435 if ((sorg == cleq->tab[j].org) && (sdest == cleq->tab[j].dest)) {
00436
00437 cleq->tab[j].org = cleq->tab[j].dest = 0 ;
00438 removed++ ;
00439 }
00440 }
00441 }
00442 }
00443
00444 qsort(&cleq->tab[0].org,
00445 cleq->numentries,
00446 sizeof(cleq_t),
00447 (qsort_fkptr_t)cleq_sort) ;
00448 if (removed < cleq->numentries) {
00449 cleq->numentries -= removed ;
00450 return removed ;
00451 } else {
00452 e_error("eliminate_redundancies: removed (%d) larger than numentries(%d)",
00453 removed, cleq->numentries) ;
00454 return -1 ;
00455 }
00456 }
00457
00458
00467
00468 static int enforce_injectiveness(cleq_tab_t * cleq)
00469 {
00470 unsigned int multmap = 0 ;
00471 size_t i,
00472 j ;
00473 label_pix minlab,
00474 dest ;
00475
00476 for (i=0 ; i<cleq->numentries ; i++) {
00477 minlab = cleq->tab[i].org ;
00478 dest = cleq->tab[i].dest ;
00479
00480 for (j=i+1 ; j<cleq->numentries ; j++)
00481 if ((dest == cleq->tab[j].dest) && (cleq->tab[j].org < minlab))
00482 minlab = cleq->tab[j].org ;
00483
00484 for (j=i ; j<cleq->numentries ; j++)
00485 if ((cleq->tab[j].org > minlab) && (dest == cleq->tab[j].dest)) {
00486 cleq->tab[j].dest = cleq->tab[j].org ;
00487 cleq->tab[j].org = minlab ;
00488 multmap ++ ;
00489 }
00490 }
00491
00492 return multmap ;
00493 }
00494
00495
00496
00504
00505 static int printcltab(cleq_tab_t *cleq)
00506 {
00507 size_t i;
00508
00509 e_comment(1, "Equivalences (%d entries):", cleq->numentries) ;
00510 for (i=0 ; i<cleq->numentries ; i++)
00511 e_comment(2, "%4d\t%4d %4d", i, cleq->tab[i].org, cleq->tab[i].dest) ;
00512 return 0 ;
00513 }
00514
00515
00524
00525 static long int relabel_image(
00526 label_image_t * l,
00527 cleq_tab_t * cl)
00528 {
00529 label_pix maxlab = 0,
00530 pix ;
00531
00532 ushort16 * cw ;
00533 label_pix * ot ;
00534 size_t nextlab = FIRST_LABEL,
00535 max_cw_index = 0 ;
00536 size_t nbpix ;
00537 size_t i,
00538 j ;
00539
00540
00541 nbpix = l->lx * l->ly ;
00542
00543
00544 cw = malloc(LUT16SIZ) ;
00545 if (cw == NULL) {
00546 e_error("relabel_image: malloc of %d bytes failed", LUT16SIZ) ;
00547 return -1 ;
00548 }
00549
00550
00551 memset(cw, NO_CORRESPONDANCE, LUT16SIZ) ;
00552
00553
00554 for (j=0 ; j<cl->numentries ; j++) {
00555 if (max_cw_index < cl->tab[j].dest) max_cw_index = cl->tab[j].dest ;
00556 cw[cl->tab[j].dest] = j ;
00557 }
00558
00559
00560 ot = calloc(CLEQSIZE, sizeof(label_pix)) ;
00561 if (ot == NULL) {
00562 e_error("relabel_image: inable to allocate ot(%d bytes)",
00563 CLEQSIZE*sizeof(label_pix)) ;
00564 free(cw) ;
00565 return -1 ;
00566 }
00567
00568
00569
00570 for (j=0 ; j<=max_cw_index ; j++) {
00571 if (cw[j] != NO_CORRESPONDANCE) pix = cl->tab[cw[j]].org ;
00572 else pix = j ;
00573 ot[pix] = pix ;
00574 if(pix > maxlab) maxlab = pix ;
00575 }
00576
00577
00578 for (i=FIRST_LABEL ; i<=maxlab ; i++) if (ot[i]) ot[i] = nextlab++ ;
00579
00580
00581 ot[(int)CNCTCMP_BLACK] = CNCTCMP_BLACK ;
00582 ot[(int)CNCTCMP_WHITE] = CNCTCMP_WHITE ;
00583
00584
00585 l->nobjtot = nextlab ;
00586 l->nrealobj = l->nobjtot - FIRST_LABEL ;
00587
00588
00589 if (!l->nrealobj)
00590 e_warning("Nothing to do (found %d real objects)", l->nrealobj) ;
00591
00592
00593 for (i=0 ; i<nbpix ; i++) {
00594 pix = l->buf[i] ;
00595 if (cw[pix] != NO_CORRESPONDANCE) l->buf[i] = ot[cl->tab[cw[pix]].org] ;
00596 else if ((pix != CNCTCMP_BLACK) && (pix != CNCTCMP_WHITE))
00597 if ((ot[pix]) && (pix != ot[pix])) {
00598 e_warning("relabel_image: cw[%d]=%d, ot[%d}=%d",
00599 pix, cw[pix], pix, ot[pix]) ;
00600 pix = ot[pix] ;
00601 l->buf[i] = pix ;
00602 }
00603 }
00604
00605
00606 free(cw) ;
00607 free(ot) ;
00608 return nextlab - FIRST_LABEL ;
00609 }
00610
00611
00612
00623
00624 static cc_stats_t *relabel_n_calc_image(
00625 label_image_t * l,
00626 cleq_tab_t * cl)
00627 {
00628 label_pix maxlab = 0,
00629 pix,
00630 maxpix ;
00631
00632 ushort16 * cw = NULL ;
00633 size_t nextlab = FIRST_LABEL,
00634 max_cw_index = 0 ;
00635 label_pix * ot ;
00636 cc_stats_t * cc = NULL,
00637 * cs = NULL ;
00638
00639 label_pix p0_1,
00640 p_10,
00641 p10,
00642 p01 ;
00643
00644 int i,
00645 j,
00646 jl ;
00647
00648
00649 if (l == NULL) {
00650 e_error("relabel_n_calc_image: NULL label_image provided") ;
00651 return NULL ;
00652 }
00653
00654
00655 cw = malloc(LUT16SIZ) ;
00656 if (cw == NULL) {
00657 e_error("relabel_n_calc_image: malloc of %d bytes failed", LUT16SIZ) ;
00658 return NULL ;
00659 }
00660
00661
00662 memset(cw, NO_CORRESPONDANCE, LUT16SIZ) ;
00663
00664
00665
00666 for (j=0 ; j<cl->numentries ; j++) {
00667 if (max_cw_index < cl->tab[j].dest) max_cw_index = cl->tab[j].dest ;
00668 cw[cl->tab[j].dest] = j ;
00669 }
00670
00671
00672 ot = calloc(CLEQSIZE, sizeof(label_pix)) ;
00673 if (ot == NULL) {
00674 e_error("relabel_n_calc_image: unable to allocate ot(%d bytes)",
00675 CLEQSIZE*sizeof(label_pix)) ;
00676 free(cw) ;
00677 return NULL ;
00678 }
00679
00680
00681
00682 for (j=0 ; j<=max_cw_index ; j++) {
00683 if (cw[j] != NO_CORRESPONDANCE) pix = cl->tab[cw[j]].org ;
00684 else pix = j ;
00685 ot[pix] = pix ;
00686 if(pix>maxlab) maxlab = pix ;
00687 }
00688
00689
00690 for(i=FIRST_LABEL ; i<=maxlab ; i++) if (ot[i]) ot[i] = nextlab++ ;
00691
00692
00693 ot[(int)CNCTCMP_BLACK] = CNCTCMP_BLACK ;
00694 ot[(int)CNCTCMP_WHITE] = CNCTCMP_WHITE ;
00695
00696
00697 l->nobjtot = nextlab ;
00698 l->nrealobj = l->nobjtot - FIRST_LABEL ;
00699
00700
00701 if (!l->nrealobj)
00702 e_warning("relabel_n_calc_image: Nothing to do (found %d real objects)",
00703 l->nrealobj) ;
00704
00705
00706 cs = new_cc_stats_table(l->nobjtot) ;
00707 if (cs == NULL) {
00708 free(ot) ;
00709 free(cw) ;
00710 return NULL ;
00711 }
00712
00713
00714 for (i=0 ; i<l->nobjtot ; i++) {
00715 cc = &cs[i] ;
00716 cc ->top_y = -1 ;
00717 cc ->rig_x = -1 ;
00718 cc ->lef_x = l->lx ;
00719 cc ->bottom_y = l->ly ;
00720 }
00721
00722
00723 maxpix = l->nobjtot ;
00724 for (j=0 ; j<l->ly ; j++) {
00725 jl = j * l->lx ;
00726 for (i=0 ; i<l->lx ; i++) {
00727 pix = l->buf[i+jl] ;
00728
00729 if (cw[pix] != NO_CORRESPONDANCE) {
00730 l->buf[i+jl] = ot[cl->tab[cw[pix]].org] ;
00731 } else {
00732 if ((pix != CNCTCMP_BLACK) && (pix != CNCTCMP_WHITE))
00733 if ((ot[pix]) && (pix != ot[pix])) {
00734 e_warning("relabel_n_calc_image: cw[%d]=%d, ot[%d}=%d",
00735 pix, cw[pix], pix, ot[pix]) ;
00736 l->buf[i+jl] = ot[pix] ;
00737 }
00738 }
00739
00740
00741 pix = l->buf[i+jl] ;
00742 if (pix >= maxpix) {
00743 e_error("inconsistent nb of real objects=%d / label=%d",
00744 l->nrealobj, pix) ;
00745 break ;
00746 }
00747 cc = &cs[pix] ;
00748
00749 cc->numpix++ ;
00750 cc->gcx += i ;
00751 cc->gcy += j ;
00752 cc->egcx += i * i ;
00753 cc->egcy += j * j ;
00754
00755
00756 get4nbs(l, i, j, &p0_1, &p_10, &p01, &p10) ;
00757 if (pix != CNCTCMP_BLACK) {
00758 if (p0_1 == CNCTCMP_BLACK) cc->wetsrf++ ;
00759 if (p_10 == CNCTCMP_BLACK) cc->wetsrf++ ;
00760 if (p01 == CNCTCMP_BLACK) cc->wetsrf++ ;
00761 if (p10 == CNCTCMP_BLACK) cc->wetsrf++ ;
00762 } else {
00763 if (p0_1 != CNCTCMP_BLACK) cc->wetsrf++ ;
00764 if (p_10 != CNCTCMP_BLACK) cc->wetsrf++ ;
00765 if (p01 != CNCTCMP_BLACK) cc->wetsrf++ ;
00766 if (p10 != CNCTCMP_BLACK) cc->wetsrf++ ;
00767 }
00768 if (j < cc->bottom_y) {
00769 cc->bottom_y = j ;
00770 cc->bottom_x = i ;
00771 }
00772 if (j > cc->top_y) {
00773 cc->top_y = j ;
00774 cc->top_x = i ;
00775 }
00776 if (i > cc->rig_x) {
00777 cc->rig_y = j ;
00778 cc->rig_x = i ;
00779 }
00780 if (i < cc->lef_x) {
00781 cc->lef_y = j ;
00782 cc->lef_x = i ;
00783 }
00784 }
00785 }
00786 free(cw) ;
00787 free(ot) ;
00788
00789
00790 for (i=0 ; i<l->nobjtot ; i++) {
00791 cc = &cs[i] ;
00792 if (cc->numpix) {
00793 if (i == CNCTCMP_WHITE) {
00794 e_warning("object=%d has size %d pixels, image badly labeled",
00795 i, cc->numpix) ;
00796 }
00797 } else {
00798 if (i != CNCTCMP_WHITE) {
00799 e_warning("object=%d has size %d pixels",
00800 i, cc->numpix) ;
00801 }
00802
00803 continue;
00804 }
00805 cc->gcx /= (double)cc->numpix ;
00806 cc->gcy /= (double)cc->numpix ;
00807 cc->egcx = sqrt(cc->egcx/(double)cc->numpix - cc->gcx) ;
00808 cc->egcy = sqrt(cc->egcy/(double)cc->numpix - cc->gcy) ;
00809 }
00810
00811 return cs ;
00812 }
00813
00814
00815
00816
00829
00830 int output_ccfeaturestable(
00831 char * fn,
00832 cc_stats_t * cs,
00833 unsigned long numobj,
00834 char * fin_name)
00835 {
00836 FILE * f = NULL ;
00837 int err = 0 ;
00838
00839 if (!numobj) e_warning("Nothing to do: numobj=%d", numobj) ;
00840 if (fn) {
00841 f = fopen(fn, "w") ;
00842 if (f == NULL) {
00843 e_error("Error opening output file %s", fn) ;
00844 return -1 ;
00845 }
00846 } else f = stdout ;
00847 err = print_ccfeaturestable(cs, numobj+2, f, fin_name) ;
00848 if (fn) fclose(f) ;
00849 return err ;
00850 }
00851
00852
00853
00866
00867 int print_ccfeaturestable(
00868 cc_stats_t * cs,
00869 unsigned long numentries,
00870 FILE * f,
00871 char * fin_name)
00872 {
00873 cc_stats_t * cc = NULL ;
00874 unsigned long j ;
00875 signed long wetsrf_chksum = 0 ;
00876
00877 if (cs == NULL) {
00878 e_error("print_ccfeaturestable: NULL cc_stats_t provided") ;
00879 return -1 ;
00880 }
00881 if (fin_name) fprintf(f,"# FILE %s.fits\n", fin_name) ;
00882 else fprintf(f,"# UNKNOWN input fits FILE\n") ;
00883 fprintf(f, "# Object=%d is Black (background)\n", CNCTCMP_BLACK) ;
00884 fprintf(f, "# Object=%d is White (unlabeled - all params should be zero)\n",
00885 CNCTCMP_WHITE);
00886 fprintf(f,
00887 "# Object size in Gravity Gravity Wet bottom top left right\n") ;
00888 fprintf(f,
00889 "# number pixels center deviation Surface\n") ;
00890 for (j=0 ; j<numentries ; j++) {
00891 cc = &cs[j] ;
00892 fprintf(f,
00893 "%8ld %8d %8.2f,%8.2f %8.2f,%8.2f %8ld %4d,%4d %4d,%4d %4d,%4d %4d,%4d\n",
00894 j,
00895 cc->numpix,
00896 (double)cc->gcx,
00897 (double)cc->gcy,
00898 (double)cc->egcx,
00899 (double)cc->egcy,
00900 cc->wetsrf,
00901 (int)cc->bottom_x,
00902 (int)cc->bottom_y,
00903 (int)cc->top_x,
00904 (int)cc->top_y,
00905 (int)cc->lef_x,
00906 (int)cc->lef_y,
00907 (int)cc->rig_x,
00908 (int)cc->rig_y) ;
00909 if (j != CNCTCMP_BLACK) wetsrf_chksum += cc->wetsrf ;
00910 }
00911 fprintf(f, "# Total %ld Objects including Black and White\n", numentries) ;
00912 fprintf(f, "# WetSurface BLACK - total is %ld\n",
00913 wetsrf_chksum - cs[(int)CNCTCMP_BLACK].wetsrf) ;
00914 return 0 ;
00915 }
00916
00917
00929
00930 static int calc_ccfeatures(
00931 label_image_t * l,
00932 cc_stats_t * cs)
00933 {
00934 label_pix pix,
00935 maxpix ;
00936
00937 label_pix p0_1,
00938 p_10,
00939 p10,
00940 p01 ;
00941 size_t i,
00942 j,
00943 jl ;
00944 cc_stats_t * cc ;
00945
00946
00947 if (l == NULL) {
00948 e_error("calc_ccfeatures: NULL label_image provided") ;
00949 return -1 ;
00950 }
00951 if (!l->nrealobj) {
00952 e_warning("calc_ccfeatures: Nothing to do (nrealobj=%d)", l->nrealobj) ;
00953 return 0 ;
00954 }
00955
00956
00957 if (cs == NULL) {
00958 e_error("calc_ccfeatures: NULL cc_stats provided") ;
00959 return -1 ;
00960 }
00961
00962
00963 maxpix = l->nobjtot ;
00964 for (j=0 ; j<l->ly ; j++) {
00965 jl = j * l->lx ;
00966 for (i=0 ; i<l->lx ; i++) {
00967 pix = l->buf[i+jl] ;
00968 if (pix >= maxpix) {
00969 e_error("inconsistent nb of objects=%d / label=%d",
00970 l->nrealobj, pix) ;
00971 break ;
00972 }
00973 cc = &cs[pix] ;
00974
00975 cc->numpix++ ;
00976 cc->gcx +=i ;
00977 cc->gcy +=j ;
00978 cc->egcx += i * i ;
00979 cc->egcy += j * j ;
00980
00981
00982 get4nbs(l, i, j, &p0_1, &p_10, &p01, &p10) ;
00983 if (pix != CNCTCMP_BLACK) {
00984 if (p0_1 == CNCTCMP_BLACK) cc->wetsrf++ ;
00985 if (p_10 == CNCTCMP_BLACK) cc->wetsrf++ ;
00986 if (p01 == CNCTCMP_BLACK) cc->wetsrf++ ;
00987 if (p10 == CNCTCMP_BLACK) cc->wetsrf++ ;
00988 } else {
00989 if (p0_1 != CNCTCMP_BLACK) cc->wetsrf++ ;
00990 if (p_10 != CNCTCMP_BLACK) cc->wetsrf++ ;
00991 if (p01 != CNCTCMP_BLACK) cc->wetsrf++ ;
00992 if (p10 != CNCTCMP_BLACK) cc->wetsrf++ ;
00993 }
00994 }
00995 }
00996
00997
00998 for (i=0 ; i<l->nobjtot ; i++) {
00999 cc = &cs[i] ;
01000 if (cc->numpix) {
01001 if (i == CNCTCMP_WHITE)
01002 e_warning("object=%d has size %d pixels, image badly labeled",
01003 i, cc->numpix) ;
01004 } else {
01005 if (i != CNCTCMP_WHITE)
01006 e_warning("calc_ccfeatures: object=%d has size %d pixels",
01007 i, cc->numpix) ;
01008
01009 continue ;
01010 }
01011 cc->gcx /= (double)cc->numpix ;
01012 cc->gcy /= (double)cc->numpix ;
01013 cc->egcx = sqrt(cc->egcx/(double)cc->numpix - cc->gcx) ;
01014 cc->egcy = sqrt(cc->egcy/(double)cc->numpix - cc->gcy) ;
01015 }
01016 return 0 ;
01017 }
01018
01019
01020
01036
01037 static label_pix search_n_classify(
01038 label_image_t * l,
01039 label_pix first_class,
01040 cleq_tab_t * cl)
01041 {
01042
01043 size_t i,
01044 j ;
01045
01046 size_t jl ;
01047
01048 int table_ok ;
01049
01050 label_pix curclass ;
01051
01052 label_pix * pdest ;
01053
01054 label_pix p0_1,
01055 p_10,
01056 p10,
01057 p01 ;
01058
01059 label_pix lb0_1,
01060 lb_10,
01061 lb00,
01062 lb10,
01063 lb01 ;
01064
01065
01066 curclass = first_class ;
01067 table_ok = TRUE ;
01068
01069
01070 for (j=0 ; j<l->ly ; j++) {
01071 jl = j * l->lx ;
01072 for (i=0 ; i<l->lx ; i++) {
01073 pdest = &l->buf[i+jl] ;
01074 if (*pdest == CNCTCMP_BLACK) continue ;
01075 if (*pdest != CNCTCMP_WHITE) {
01076 if (curclass <= *pdest) curclass = *pdest + 1 ;
01077 if (table_ok) table_ok = add_equiv_cl(cl, *pdest, *pdest) ;
01078 continue ;
01079 }
01080 get4nbs(l, i, j, &p0_1, &p_10, &p10, &p01) ;
01081 lb_10 = lb0_1 = lb10 = lb01 = 0 ;
01082
01083
01084 switch(p_10) {
01085 case CNCTCMP_BLACK:
01086
01087 break ;
01088
01089 case CNCTCMP_WHITE:
01090
01091 e_error("labeling order error at p_10(%d,%d)=%d", i, j, p_10) ;
01092 return 0 ;
01093
01094 default:
01095 lb_10 = p_10 ;
01096 }
01097
01098 switch(p0_1) {
01099 case CNCTCMP_BLACK:
01100
01101 break ;
01102
01103 case CNCTCMP_WHITE:
01104
01105 e_error("labeling order error at p0_1(%d,%d)", i, j) ;
01106 return 0 ;
01107
01108 default:
01109 if (p0_1 != lb_10) lb0_1 = p0_1 ;
01110 }
01111
01112 switch(p01) {
01113
01114 case CNCTCMP_BLACK:
01115 case CNCTCMP_WHITE:
01116 break ;
01117 default:
01118 if ((p01 != lb_10) && (p01 !=lb0_1)) lb01 = p01 ;
01119 }
01120
01121 switch(p10){
01122
01123 case CNCTCMP_BLACK:
01124 case CNCTCMP_WHITE:
01125 break ;
01126 default:
01127 if ((p10 != lb_10) && (p10 !=lb0_1) && (p10 != lb01)) lb10 = p10 ;
01128 }
01129
01130
01131 if (*pdest == CNCTCMP_WHITE) {
01132 if (lb_10) *pdest = lb_10 ;
01133 else if (lb0_1) *pdest = lb0_1 ;
01134 else if (lb01) *pdest = lb01 ;
01135 else if (lb10) *pdest = lb10 ;
01136 else {
01137 if (curclass < MAX_LABEL) {
01138
01139 if (curclass != CNCTCMP_WHITE) *pdest = curclass ;
01140 else e_error("labeling order error at p00(%d,%d)=%d lb_10=%d, lb0_1=%d, lb01=%d, lb10=%d, cc=%d",
01141 i,j,*pdest,lb_10,lb0_1,lb01,curclass) ;
01142 if (table_ok) table_ok = add_equiv_cl(cl,
01143 curclass,
01144 curclass) ;
01145 curclass++ ;
01146 } else {
01147 table_ok = FALSE ;
01148 break ;
01149 }
01150 }
01151 }
01152
01153 if ((*pdest != lb_10) && (*pdest !=lb0_1) &&
01154 (*pdest != lb01) && (*pdest !=lb10)) lb00 = *pdest ;
01155 else lb00 = CNCTCMP_BLACK ;
01156
01157
01158 if (table_ok) table_ok = add_equiv_cl(cl, lb00, lb0_1) ;
01159 if (table_ok) table_ok = add_equiv_cl(cl, lb00, lb01) ;
01160 if (table_ok) table_ok = add_equiv_cl(cl, lb00, lb10) ;
01161 if (table_ok) table_ok = add_equiv_cl(cl, lb00, lb_10) ;
01162 if (table_ok) table_ok = add_equiv_cl(cl, lb_10, lb0_1) ;
01163 if (table_ok) table_ok = add_equiv_cl(cl, lb_10, lb01) ;
01164 if (table_ok) table_ok = add_equiv_cl(cl, lb_10, lb10) ;
01165 if (table_ok) table_ok = add_equiv_cl(cl, lb0_1, lb01) ;
01166 if (table_ok) table_ok = add_equiv_cl(cl, lb0_1, lb10) ;
01167 if (table_ok) table_ok = add_equiv_cl(cl, lb01, lb10) ;
01168
01169
01170 if (!table_ok) break;
01171 }
01172
01173
01174 if (!table_ok) break ;
01175 }
01176
01177
01178 return curclass - first_class ;
01179 }
01180
01181
01182
01209
01210 cc_stats_t * labelize_binary(label_image_t * label_image)
01211 {
01212 cc_stats_t * cs = NULL ;
01213 cleq_tab_t * cleq ;
01214 unsigned long new_classes = 0 ;
01215 label_pix numobj = 0 ;
01216 int Cycle = 1,
01217 redundancies,
01218 equival,
01219 reductions ;
01220 int Feature_calc_Done = FALSE ;
01221 int err = 0 ;
01222
01223 if (label_image == NULL) {
01224 e_error("NULL label_image: aborting labelize_binary") ;
01225 return NULL ;
01226 }
01227 if (label_image->buf == NULL) {
01228 e_error("NULL label_image buffer: aborting labelize_binary") ;
01229 return NULL ;
01230 }
01231
01232 cleq = calloc(1, sizeof(cleq_tab_t)) ;
01233 if (cleq == NULL) {
01234 e_error("labelize_binary: unable to allocate cleq (%d bytes)",
01235 sizeof(cleq_tab_t)) ;
01236 return NULL ;
01237 }
01238
01239 cleq->tab = calloc(CLEQSIZE, sizeof(cleq_t)) ;
01240 if (cleq->tab == NULL) {
01241 e_error("labelize_binary: unable to allocate cleq->tab (%d bytes)",
01242 CLEQSIZE*sizeof(cleq_t)) ;
01243 free(cleq) ;
01244 return NULL ;
01245 }
01246
01247 do {
01248 new_classes = search_n_classify(label_image, FIRST_LABEL, cleq ) ;
01249 do {
01250 equival = resolve_transitivities(cleq) ;
01251 redundancies = eliminate_redundancies(cleq) ;
01252 if (redundancies < 0) {
01253 err = redundancies ;
01254 printcltab(cleq) ;
01255 break ;
01256 }
01257
01258 reductions = enforce_injectiveness(cleq) ;
01259 } while ((reductions || redundancies) && !err) ;
01260 if (new_classes >= MAX_CLASSNUM) {
01261 e_comment(2, "pass %d: merging & reordering equivalent clusters", 2) ;
01262 numobj = relabel_image(label_image, cleq) ;
01263 memset(cleq->tab, 0, sizeof(cleq_t)*CLEQSIZE) ;
01264 cleq->numentries = 0 ;
01265 } else {
01266 cs = relabel_n_calc_image(label_image, cleq) ;
01267 numobj = label_image->nrealobj ;
01268 Feature_calc_Done = TRUE ;
01269 }
01270 if (new_classes >= MAX_CLASSNUM) e_comment(2, "cycle %d:", ++Cycle) ;
01271 } while (new_classes && new_classes >= MAX_CLASSNUM && !err && (Cycle<20)) ;
01272
01273 if (!Feature_calc_Done) {
01274 cs = new_cc_stats_table(numobj) ;
01275 err = calc_ccfeatures(label_image, cs) ;
01276 }
01277 label_image->nrealobj = numobj ;
01278 free(cleq->tab) ;
01279 free(cleq) ;
01280
01281 if (!err) return cs ;
01282 else return NULL ;
01283 }
01284
01285
01296
01297 int dump_labelimage(
01298 label_image_t * to_dump,
01299 char * filename)
01300 {
01301 FILE * out ;
01302 fits_header * fh ;
01303 size_t to_write,
01304 written ;
01305 char cval[80] ;
01306
01307 if (filename == NULL) {
01308 e_error("no file name given: aborting dumping of label image") ;
01309 return -1 ;
01310 }
01311 if (to_dump == NULL) {
01312 e_error("NULL input label image: aborting dump") ;
01313 return -1 ;
01314 }
01315
01316 e_comment(2, "saved %8d real objects, total labels %d",
01317 to_dump->nrealobj, to_dump->nobjtot) ;
01318 if (to_dump->nrealobj < 1)
01319 e_warning("Warning! No objects found in label image file %s", filename) ;
01320
01321 fh = fits_header_default() ;
01322 sprintf(cval, "%d", BPP_16_SIGNED) ;
01323 fits_header_add(fh, "BITPIX", cval, "bits per pixel", NULL) ;
01324 sprintf(cval, "2") ;
01325 fits_header_add(fh, "NAXIS", cval, "single image", NULL) ;
01326 sprintf(cval, "%d", to_dump->lx) ;
01327 fits_header_add(fh, "NAXIS1", cval, "x axis", NULL) ;
01328 sprintf(cval, "%d", to_dump->ly) ;
01329 fits_header_add(fh, "NAXIS2", cval, "x axis", NULL) ;
01330
01331
01332 sprintf(cval, "%d", to_dump->nobjtot) ;
01333 fits_header_add(fh, "LABELS", cval, "total objects", NULL) ;
01334
01335 sprintf(cval, "%d", to_dump->nrealobj) ;
01336 fits_header_add(fh, "REALOBJS", cval, "real objects", NULL) ;
01337
01338 if ((out = fopen(filename, "a")) == NULL) {
01339 e_error("cannot append dump file [%s]", filename) ;
01340 return -1 ;
01341 }
01342 fits_dump_header(fh, out) ;
01343 fits_header_destroy(fh) ;
01344
01345 to_write = (size_t)(to_dump->lx * to_dump->ly) ;
01346 written = fwrite(to_dump->buf,BYTESPERPIXEL(BPP_16_SIGNED), to_write, out) ;
01347 fclose(out) ;
01348
01349 if (to_write != written) {
01350 e_error("%ld bytes written to disk instead of %ld", written, to_write) ;
01351 e_error("cannot write data to disk (disk full?): aborting") ;
01352 return -1 ;
01353 }
01354
01355 stuff_with_zeros(filename, to_dump->lx * to_dump->ly) ;
01356 return 0 ;
01357 }
01358
01359
01360
01369
01370 static int get_fits_info(
01371 char * filename,
01372 size_t * lx,
01373 size_t * ly,
01374 int * ptype,
01375 size_t * headersize)
01376 {
01377 char * cval = NULL ;
01378
01379 if (lx == NULL) {
01380 e_error("get_fits_info: lx = %d", lx) ;
01381 return -1 ;
01382 }
01383 if (ly == NULL) {
01384 e_error("get_fits_info: ly = %d", ly) ;
01385 return -1 ;
01386 }
01387 if (ptype == NULL) {
01388 e_error("get_fits_info: ptype NULL") ;
01389 return -1 ;
01390 }
01391 if (headersize == NULL){
01392 e_error("get_fits_info: headersize NULL") ;
01393 return -1 ;
01394 }
01395
01396
01397 if (eclipse_is_fits_file(filename) != 1) {
01398 e_error("file [%s] is not a pixel map", filename) ;
01399 return -1 ;
01400 }
01401
01402
01403 cval = query_fits_header(filename, "NAXIS1") ;
01404 *lx = atoi(cval) ;
01405 if ((*lx < 1) || (*lx > MAX_COLUMN_NUMBER)) {
01406 e_error("cannot process pixel map with NAXIS1 = %d", lx) ;
01407 return -1 ;
01408 }
01409
01410
01411 cval = query_fits_header(filename, "NAXIS2") ;
01412 *ly = atoi(cval) ;
01413 if ((*ly < 1) || (*ly > MAX_LINE_NUMBER)) {
01414 e_error("cannot process pixel map with NAXIS2 = %d", ly) ;
01415 return -1 ;
01416 }
01417
01418
01419 cval = query_fits_header(filename, "BITPIX") ;
01420 *ptype = atoi(cval) ;
01421
01422 return 0 ;
01423 }
01424
01425
01426
01437
01438 label_image_t * load_label_image_from_pixelmap(char * filename)
01439 {
01440 label_image_t * lab ;
01441 pixel_map * map ;
01442 int i ;
01443
01444 map = load_pixelmap(filename) ;
01445 if (map == NULL) {
01446 e_error("cannot load pixelmap [%s]", filename) ;
01447 return NULL ;
01448 }
01449 lab = new_label_image(map->lx, map->ly) ;
01450 if (lab == NULL) {
01451 e_error("label_image allocation failure: aborting load") ;
01452 return NULL ;
01453 }
01454 for (i=0 ; i<map->nbpix ; i++) {
01455 lab->buf[i] = (label_pix)map->data[i] ;
01456 }
01457 destroy_pixelmap(map) ;
01458 return lab ;
01459 }
01460
01461
01473
01474 label_image_t * pixelmap_2_label_image(pixel_map *pm)
01475 {
01476 label_image_t * lab ;
01477 int i ;
01478
01479 if (pm == NULL) return NULL ;
01480
01481 lab = new_label_image(pm->lx, pm->ly) ;
01482 for (i=0 ; i<pm->nbpix ; i++) {
01483 lab->buf[i] = (label_pix)pm->data[i] ;
01484 }
01485 return lab ;
01486 }
01487
01488
01496
01497 label_image_t * load_label_image(char * filename)
01498 {
01499 FILE * in = NULL ;
01500 label_image_t * label_image = NULL;
01501 size_t lx,
01502 ly ;
01503 int ptype,
01504 err = 0 ;
01505 size_t headersize =0 ;
01506
01507 err = get_fits_info(filename, &lx, &ly, &ptype, &headersize) ;
01508 if (err) return NULL ;
01509
01510 label_image = new_label_image(lx, ly) ;
01511 if (label_image == NULL) {
01512 e_error("label_image allocation failure: aborting load") ;
01513 return NULL ;
01514 }
01515
01516 if ((in = fopen(filename, "r")) == NULL) {
01517 e_error("cannot access file %s any more: aborting load", filename) ;
01518 free_label_image(label_image) ;
01519 return NULL ;
01520 }
01521
01522 fseek(in, headersize, SEEK_SET) ;
01523 fread(label_image->buf, lx*ly, BYTESPERPIXEL(BPP_16_SIGNED), in) ;
01524 fclose(in) ;
01525
01526 return label_image ;
01527 }
01528
01529
01547
01548 int get_extreme_obj_coor(
01549 cc_stats_t * cs,
01550 int nobjs,
01551 int * selected_objlist,
01552 int nselectedobjs,
01553 int * left,
01554 int * right,
01555 int * top,
01556 int * bot)
01557 {
01558 int err = 0 ;
01559 int i,
01560 l,
01561 r,
01562 b,
01563 t ;
01564
01565 if (cs == NULL) {
01566 e_error("get_extreme_obj_coor: objects stats NULL") ;
01567 return -1 ;
01568 }
01569 if (nobjs <= 0) {
01570 e_error("get_extreme_obj_coor: nobjs %d",nobjs) ;
01571 return -1 ;
01572 }
01573 if (selected_objlist == NULL) {
01574 l = cs[FIRST_LABEL].lef_x ;
01575 r = cs[FIRST_LABEL].rig_x ;
01576 t = cs[FIRST_LABEL].top_y ;
01577 b = cs[FIRST_LABEL].bottom_y ;
01578 for (i=FIRST_LABEL ; i<nobjs ; i++) {
01579 if (cs[i].lef_x < l) l = cs[i].lef_x ;
01580 if (cs[i].rig_x > r) r = cs[i].rig_x ;
01581 if (cs[i].bottom_y < b) b = cs[i].bottom_y ;
01582 if (cs[i].top_y > t) t = cs[i].top_y ;
01583 }
01584 } else {
01585 if (!nselectedobjs) return -1 ;
01586 i = selected_objlist[nselectedobjs-1] ;
01587 l = cs[i].lef_x ;
01588 r = cs[i].rig_x ;
01589 t = cs[i].top_y ;
01590 b = cs[i].bottom_y ;
01591 while (--nselectedobjs) {
01592 i = selected_objlist[nselectedobjs] ;
01593 if (cs[i].lef_x < l) l = cs[i].lef_x ;
01594 if (cs[i].rig_x > r) r = cs[i].rig_x ;
01595 if (cs[i].bottom_y < b) t = cs[i].bottom_y ;
01596 if (cs[i].top_y > t) t = cs[i].top_y ;
01597 }
01598 }
01599
01600 *left = l ;
01601 *right = r ;
01602 *bot = b ;
01603 *top = t ;
01604
01605
01606
01607 if (debug_active())
01608 e_comment(0,
01609 "get_extreme_obj_coor: left=%d, right=%d, top=%d, bottom=%d", l, r, t, b) ;
01610
01611
01612 return err ;
01613 }
01614
01615
01625
01626 int select_objs(
01627 label_image_t * l,
01628 int * objlist,
01629 int nobjs)
01630 {
01631 int err = 0 ;
01632 long int i,
01633 nbpix ;
01634 int * lut = NULL ;
01635
01636 if (l == NULL) {
01637 e_error("select_objs : label image NULL") ;
01638 return -1 ;
01639 }
01640 if (objlist == NULL) {
01641 e_error("select_objs : objlist NULL") ;
01642 return -1 ;
01643 }
01644 if (nobjs > l->nrealobj) {
01645 e_error("select_objs : objects know %d objects to select %d",
01646 l->nrealobj, nobjs) ;
01647 return -1 ;
01648 }
01649 nbpix = l->lx * l->ly ;
01650 if (!nbpix) {
01651 e_error("select_objs : nbpix %d", nbpix) ;
01652 return -1 ;
01653 }
01654 lut = calloc(l->nobjtot, sizeof(int)) ;
01655 if (lut == NULL) {
01656 e_error("select_objs : failed allocating lut for %d objects", nobjs) ;
01657 return -1 ;
01658 }
01659 for (i=0 ; i<nobjs ; i++) lut[objlist[i]] = objlist[i] ;
01660
01661 e_comment(0, "select_objs : deleting %d objects", l->nrealobj - nobjs) ;
01662 l->nrealobj = nobjs ;
01663 for (i=0;i<nbpix;i++) l->buf[i] = lut[l->buf[i]] ;
01664 free(lut) ;
01665
01666 return err ;
01667 }
01668
01669
01670
01682
01683 int mask_objs(
01684 OneImage * in,
01685 OneImage * out,
01686 label_image_t * l,
01687 int * objlist,
01688 int nobjs,
01689 pixelvalue maskval)
01690 {
01691 int err = 0 ;
01692 int i,
01693 nbpix ;
01694 int * lut ;
01695
01696 if (l == NULL) {
01697 e_error("mask_objs: label image NULL") ;
01698 return -1 ;
01699 }
01700 if (objlist == NULL) {
01701 e_error("mask_objs: objlist NULL") ;
01702 return -1 ;
01703 }
01704 if (nobjs > l->nrealobj) {
01705 e_error("mask_objs: objects known %d, objects to select %d",
01706 l->nrealobj, nobjs) ;
01707 return -1 ;
01708 }
01709 nbpix = l->lx * l->ly ;
01710 if (nbpix <= 0) {
01711 e_error("select_objs : nbpix %d", nbpix) ;
01712 return -1 ;
01713 }
01714 lut = calloc(l->nobjtot, sizeof(int)) ;
01715 if (lut == NULL) {
01716 e_error("select_objs : failed allocating lut for %d objects", nobjs) ;
01717 return -1 ;
01718 }
01719 for (i=0 ; i<nobjs ; i++) lut[objlist[i]] = objlist[i] ;
01720
01721 for (i=0 ; i<nbpix ; i++) {
01722 if (lut[l->buf[i]]) out->data[i] = in->data[i] ;
01723 else out->data[i] = maskval ;
01724 }
01725 free(lut) ;
01726
01727 return err ;
01728 }
01729
01730
01731
01732
01733