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 }