my_cnct_comp.c

00001 /*-------------------------------------------------------------------------
00002    
00003    File name    :   cnct_comp.c
00004    Author       :   Thomas Rogon
00005    Created on   :   13 April 1999
00006    Description  :   connected components handling
00007 
00008  *--------------------------------------------------------------------------*/
00009 /*
00010     $Id: cnct_comp.c,v 1.32 2000/10/11 14:09:05 ndevilla Exp $
00011     $Author: ndevilla $
00012     $Date: 2000/10/11 14:09:05 $
00013     $Revision: 1.32 $
00014 */
00015 
00016 /*---------------------------------------------------------------------------
00017                                 Includes
00018  ---------------------------------------------------------------------------*/
00019 #include "eclipse.h"
00020 #include "cnct_comp.h"
00021 
00022 
00023 /*---------------------------------------------------------------------------
00024                                 Defines
00025  ---------------------------------------------------------------------------*/
00026 #define CLEQSIZE        65534 /* limit class equivalence table size */
00027 #define MAXUI16         65535
00028 
00029 
00030 #define LUT16SIZ    ((MAXUI16+1)*sizeof(ushort16))
00031 #define NO_CORRESPONDANCE (0xFFFF)
00032 
00033 
00034 /*---------------------------------------------------------------------------
00035                                 New types
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                             Function prototypes
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                             Function codes
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     /* Check the requested size */
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     /* Set up nbpix */
00212     nbpix = lx * ly ;
00213 
00214     /* Memory allocation */
00215     lab = (label_image_t*)calloc(1, sizeof(label_image_t)) ;
00216     
00217     /* Fill the fields */
00218     lab->lx = lx ;
00219     lab->ly = ly ;
00220     lab->nrealobj = 0 ;
00221     lab->nobjtot = 0 ;
00222 
00223     /* Pixels allocation */
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     /* Return */
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     /* Check the requested size validity */
00252     if (!numobj) e_warning("new_cc_stats: nothing to do (numobj=%d)", numobj) ;
00253 
00254     /* Memory allocation */
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     /* Return */
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     /* zero entries = empty */
00326     while ((i < CLEQSIZE) && cl->tab[i].org) { 
00327         /*entry already in table */
00328         if ((cl->tab[i].org == sorg) && (cl->tab[i].dest == sdest)) return TRUE ; 
00329         else i++ ;
00330     }
00331     
00332     /* add entry */
00333     if (i < CLEQSIZE) { 
00334         /* consistency check */
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     /* table saturated */
00348     return FALSE ; 
00349 }
00350 
00351 
00352 /*---------------------------------------------------------------------------
00353  * For qsort(...) 
00354  * Notice: valid entries are nonzero 
00355  *--------------------------------------------------------------------------*/
00356 static int cleq_sort(
00357         const cleq_t    *   c1, 
00358         const cleq_t    *   c2)
00359 { 
00360     /* both valid */
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    /* void c1 always yields to valid c2 */
00368    if (c2->org) return 1 ;
00369    /* void c2 always yields to valid c1 */
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                     /* clear entry */
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     /* correspondance table  16 -> 16 bits */
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     /* 16 bit lookup of 16 bit */
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     /* mark unused entries */
00551     memset(cw, NO_CORRESPONDANCE, LUT16SIZ) ;
00552     
00553     /* Find out the index max and fill cw: cw[dest] = index */
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     /* Allocate 'output table' : ot[dest] = org  */
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     /* Find out maxlab and set ot : ot[i]=i for each i not in destination */
00569     /* ot[something present in destination] = 0 */
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     /* ot[2]=2,ot[3]=0,ot[4]=4,ot[5]=5 --> ot[2]=2,ot[3]=0,ot[4]=3, ot[5]=4 */
00578     for (i=FIRST_LABEL ; i<=maxlab ; i++) if (ot[i]) ot[i] = nextlab++ ;
00579     
00580     /* Special cases */
00581     ot[(int)CNCTCMP_BLACK] = CNCTCMP_BLACK ;
00582     ot[(int)CNCTCMP_WHITE] = CNCTCMP_WHITE ;
00583 
00584     /* Update label_image object */
00585     l->nobjtot = nextlab ;
00586     l->nrealobj = l->nobjtot - FIRST_LABEL ;
00587 
00588     /* No object left in the image */
00589     if (!l->nrealobj)
00590         e_warning("Nothing to do (found %d real objects)", l->nrealobj) ;
00591 
00592     /* Modify the input label image if necessary     */
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     /* Free and return */
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     /* correspondance table  16 -> 16 bits */
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     /* 4 neighbour pixels */
00639     label_pix       p0_1,
00640                     p_10,   
00641                     p10,
00642                     p01 ;
00643 
00644     int             i,
00645                     j,
00646                     jl ;
00647     
00648     /* Test input label image */
00649     if (l == NULL) {
00650         e_error("relabel_n_calc_image: NULL label_image provided") ;
00651         return NULL ;
00652     }
00653 
00654     /* 16 bit lookup of 16 bit */
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     /* Mark unused entries */
00662     memset(cw, NO_CORRESPONDANCE, LUT16SIZ) ; 
00663     
00664     
00665     /* Find out the index max and fill cw: cw[dest] = index */
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     /* Allocate 'output table' : ot[dest] = org  */
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     /* Find out maxlab and set ot : ot[i]=i for each i not in destination */
00681     /* ot[something present in destination] = 0 */
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     /* ot[2]=2,ot[3]=0,ot[4]=4,ot[5]=5 --> ot[2]=2,ot[3]=0,ot[4]=3, ot[5]=4 */
00690     for(i=FIRST_LABEL ; i<=maxlab ; i++) if (ot[i]) ot[i] = nextlab++ ;
00691     
00692     /* Special cases */
00693     ot[(int)CNCTCMP_BLACK] = CNCTCMP_BLACK ;
00694     ot[(int)CNCTCMP_WHITE] = CNCTCMP_WHITE ;
00695 
00696     /* Update label_image object */
00697     l->nobjtot = nextlab ;
00698     l->nrealobj = l->nobjtot - FIRST_LABEL ;
00699 
00700     /* No object left in the image */
00701     if (!l->nrealobj)
00702         e_warning("relabel_n_calc_image: Nothing to do (found %d real objects)",
00703                     l->nrealobj) ;
00704 
00705     /* New statistic table */
00706     cs = new_cc_stats_table(l->nobjtot) ;
00707     if (cs == NULL) {
00708         free(ot) ;
00709         free(cw) ;
00710         return NULL ;
00711     }
00712     
00713     /* Initalise features */
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     /* Computing features for l->nrealobj real objects */
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             /* Modify the input label image if necessary     */
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             /* Calc features */
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             /* calculate wet surface by considering neighbour configurations*/          
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     /* Sum it up */
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             /* don't divide by size 0*/
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     /* neighbour pixels */
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     /* Test on input label image */
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     /* Test on stat object */
00957     if (cs == NULL) {
00958         e_error("calc_ccfeatures: NULL cc_stats provided") ;
00959         return -1 ;
00960     }
00961     
00962     /* Computing features */
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             /* Calculate wet surface by considering neighbour configurations*/          
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     /* Sum it up */
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             /* don't divide by size 0*/
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     /* row and col   */
01043     size_t          i,
01044                     j ; 
01045     /* row*col combinations for faster indexing in image */
01046     size_t          jl ; 
01047     /* flag signalling class equiavlence table not saturated */
01048     int             table_ok ; 
01049     /* intermediate class counter */
01050     label_pix       curclass ;
01051     /* pointer to destinatation pixel */
01052     label_pix   *   pdest ; 
01053     /* pixel neighbourhood definition */
01054     label_pix       p0_1, 
01055                     p_10,
01056                     p10,
01057                     p01 ;
01058     /* labels assigned within neighbourhood */
01059     label_pix       lb0_1,
01060                     lb_10, 
01061                     lb00,
01062                     lb10,
01063                     lb01 ;
01064 
01065     /* Initialise */
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             /* check defined neighbours*/
01084             switch(p_10) {
01085                 case CNCTCMP_BLACK:
01086                 /* skip */
01087                 break ;
01088                 
01089                 case CNCTCMP_WHITE: 
01090                 /* preceding point should have been labelled */
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                 /*skip*/
01101                 break ;
01102          
01103                 case CNCTCMP_WHITE: 
01104                 /*preceding point should have been labelled*/
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                 /*unlabelled*/
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                 /*unlabelled*/
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             /* label point */
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                         /* avoid using rsvd label*/
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             /* add equivalences, avoid redundant relations*/
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             /* terminate if table saturated! */
01170             if (!table_ok) break; 
01171         }           
01172         
01173         /* terminate if table saturated! */
01174         if (!table_ok) break ;
01175     }
01176 
01177     /* report nb of new classes found */
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     /* Add up number of labels */
01332     sprintf(cval, "%d", to_dump->nobjtot) ;
01333     fits_header_add(fh, "LABELS", cval, "total objects", NULL) ;
01334     /* Add number of real objects */ 
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     /* Test file existence and FITSability  */
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     /* Get First axis size and check it */
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     /* Get 2nd axis size and check it   */
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     /* Test data type */
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     /* DEBUG CASE */
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 /*--------------------------------------------------------------------------*/

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