00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include "mrilib.h"
00008 
00009 /*----------------------------------------------------------------
00010    Find connected components in a graph.
00011    Inputs:
00012      int    npt  = # of points
00013      int ** gmat = gmat[i] is the connectivity list of pt # i
00014                    gmat[i][0] (=nc) is the number of connections
00015                    (if 0, can have gmat[i] = NULL to save space)
00016                    gmat[i][1..nc] is the index of connected points
00017                    (if gmat[i][j] < 0, this pt is skipped)
00018    Output
00019      int * ncom   = *ncom is the number of components found
00020      int *** cmat = *cmat[j] is the list of points in component # j
00021                     *cmat[j][0] (=nc) is the number of points
00022                     *cmat[j][1..nc] is the list of points
00023 
00024    If *ncom is returned as 0, something bad happened.  Otherwise,
00025    *cmat and *cmat[j] are malloc()-ed and should be free()-ed
00026    someday.
00027 ------------------------------------------------------------------*/
00028 
00029 void GRAPH_find_components( int npt, int ** gmat, int * ncom, int *** cmat )
00030 {
00031    int ** cm ;  /* will be *cmat */
00032    int ncm ;    /* will be *ncom */
00033    int cmall , cmuse ;
00034    byte * used ; /* used up point list */
00035    int ijk , ijk_last , ncon , ii,kk,pkk,nkk,pii , *tcm ;
00036 
00037    if( ncom == NULL ) return ;                             /* error */
00038    *ncom = 0 ;                                             /* default return */
00039    if( npt <= 0 || gmat == NULL || cmat == NULL ) return ; /* error */
00040 
00041    cm  = (int **) malloc( sizeof(int *) * 1 ) ;            /* initialize */
00042    ncm = 0 ;                                               /* # of compon */
00043 
00044    used = (byte *) calloc( npt , sizeof(byte) ) ;          /* used up? */
00045 
00046    /*--- scan through array,
00047          find nonzero point,
00048          build a component about it, start scan over again ---*/
00049 
00050    ijk_last = 0 ;
00051    do {
00052       for( ijk=ijk_last ; ijk < npt ; ijk++ )
00053          if( gmat[ijk] != NULL && gmat[ijk][0] > 0 && !used[ijk] ) break ;
00054 
00055       if( ijk == npt ) break ;  /* didn't find any! */
00056 
00057       ijk_last = ijk+1 ;        /* start here next time */
00058 
00059       /* make a new component */
00060 
00061       used[ijk] = 1 ;                           /* mark this point as used */
00062       ncon = gmat[ijk][0] ;
00063 
00064       ncm++ ;                                                /* # compon   */
00065       cm = (int **) realloc( cm , sizeof(int *) * ncm ) ;    /* add compon  */
00066       cm[ncm-1] = (int *) malloc( sizeof(int) * (ncon+8) ) ; /* compon array */
00067       cmuse = 1 ;                                            /* has 1 point   */
00068       cm[ncm-1][1] = ijk ;                                   /* add that point */
00069       cmall = (ncon+8) ;                                     /* space allocated */
00070 
00071       /*--
00072         for each point now in component:
00073            add all points connected to it to the component,
00074              that aren't already used up
00075            continue until end of component is reached
00076              (note that component is expanding as we progress)
00077       --*/
00078 
00079       for( kk=1 ; kk <= cmuse ; kk++ ){
00080          pkk = cm[ncm-1][kk] ;                 /* kk-th point in component */
00081          nkk = (gmat[pkk] == NULL) ? 0         /* # pts attached to it */
00082                                    : gmat[pkk][0] ;
00083 
00084          for( ii=1 ; ii <= nkk ; ii++ ){
00085             pii = gmat[pkk][ii] ;              /* ii-th point attached to pkk */
00086             if( pii >= 0 && !used[pii] ){      /* pii not used yet? */
00087                                                /* then add pii to component */
00088 
00089                if( ++cmuse >= cmall ){         /* expand component if needed */
00090                   cmall     = cmuse+32 ;
00091                   cm[ncm-1] = (int *) realloc( cm[ncm-1], sizeof(int)*cmall ) ;
00092                }
00093 
00094                cm[ncm-1][cmuse] = pii ;        /* add pii to component */
00095                used[pii] = 1 ;                 /* mark pii as used */
00096             }
00097          }
00098       }                                        /* this component is done */
00099 
00100       cm[ncm-1][0] = cmuse ;                   /* store size of component */
00101 
00102       if( cmall > cmuse+1 )                    /* trim component array */
00103          cm[ncm-1] = (int *) realloc( cm[ncm-1], sizeof(int)*(cmuse+1) ) ;
00104 
00105    } while( 1 ) ;  /* end of loop over components */
00106 
00107    /* prepare for output */
00108 
00109    free(used) ;                                /* toss the trash */
00110 
00111    if( ncm == 0 ){ free(cm) ; cm = NULL ; }    /* no components!? */
00112 
00113    /* sort components by size */
00114 
00115    for( kk=0 ; kk < ncm ; kk++ ){
00116       for( ijk=0,ii=1 ; ii < ncm ; ii++ ){
00117          if( cm[ii-1][0] < cm[ii][0] ){
00118             tcm = cm[ii-1] ; cm[ii-1] = cm[ii] ; cm[ii] = tcm ; ijk++ ;
00119          }
00120       }
00121       if( !ijk ) break ;
00122    }
00123 
00124    *ncom = ncm ; *cmat = cm ; return ;
00125 }