1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8
9 /*----------------------------------------------------------------
10 Find connected components in a graph.
11 Inputs:
12 int npt = # of points
13 int ** gmat = gmat[i] is the connectivity list of pt # i
14 gmat[i][0] (=nc) is the number of connections
15 (if 0, can have gmat[i] = NULL to save space)
16 gmat[i][1..nc] is the index of connected points
17 (if gmat[i][j] < 0, this pt is skipped)
18 Output
19 int * ncom = *ncom is the number of components found
20 int *** cmat = *cmat[j] is the list of points in component # j
21 *cmat[j][0] (=nc) is the number of points
22 *cmat[j][1..nc] is the list of points
23
24 If *ncom is returned as 0, something bad happened. Otherwise,
25 *cmat and *cmat[j] are malloc()-ed and should be free()-ed
26 someday.
27 ------------------------------------------------------------------*/
28
GRAPH_find_components(int npt,int ** gmat,int * ncom,int *** cmat)29 void GRAPH_find_components( int npt, int ** gmat, int * ncom, int *** cmat )
30 {
31 int ** cm ; /* will be *cmat */
32 int ncm ; /* will be *ncom */
33 int cmall , cmuse ;
34 byte * used ; /* used up point list */
35 int ijk , ijk_last , ncon , ii,kk,pkk,nkk,pii , *tcm ;
36
37 if( ncom == NULL ) return ; /* error */
38 *ncom = 0 ; /* default return */
39 if( npt <= 0 || gmat == NULL || cmat == NULL ) return ; /* error */
40
41 cm = (int **) malloc( sizeof(int *) * 1 ) ; /* initialize */
42 ncm = 0 ; /* # of compon */
43
44 used = (byte *) calloc( npt , sizeof(byte) ) ; /* used up? */
45
46 /*--- scan through array,
47 find nonzero point,
48 build a component about it, start scan over again ---*/
49
50 ijk_last = 0 ;
51 do {
52 for( ijk=ijk_last ; ijk < npt ; ijk++ )
53 if( gmat[ijk] != NULL && gmat[ijk][0] > 0 && !used[ijk] ) break ;
54
55 if( ijk == npt ) break ; /* didn't find any! */
56
57 ijk_last = ijk+1 ; /* start here next time */
58
59 /* make a new component */
60
61 used[ijk] = 1 ; /* mark this point as used */
62 ncon = gmat[ijk][0] ;
63
64 ncm++ ; /* # compon */
65 cm = (int **) realloc( cm , sizeof(int *) * ncm ) ; /* add compon */
66 cm[ncm-1] = (int *) malloc( sizeof(int) * (ncon+8) ) ; /* compon array */
67 cmuse = 1 ; /* has 1 point */
68 cm[ncm-1][1] = ijk ; /* add that point */
69 cmall = (ncon+8) ; /* space allocated */
70
71 /*--
72 for each point now in component:
73 add all points connected to it to the component,
74 that aren't already used up
75 continue until end of component is reached
76 (note that component is expanding as we progress)
77 --*/
78
79 for( kk=1 ; kk <= cmuse ; kk++ ){
80 pkk = cm[ncm-1][kk] ; /* kk-th point in component */
81 nkk = (gmat[pkk] == NULL) ? 0 /* # pts attached to it */
82 : gmat[pkk][0] ;
83
84 for( ii=1 ; ii <= nkk ; ii++ ){
85 pii = gmat[pkk][ii] ; /* ii-th point attached to pkk */
86 if( pii >= 0 && !used[pii] ){ /* pii not used yet? */
87 /* then add pii to component */
88
89 if( ++cmuse >= cmall ){ /* expand component if needed */
90 cmall = cmuse+32 ;
91 cm[ncm-1] = (int *) realloc( cm[ncm-1], sizeof(int)*cmall ) ;
92 }
93
94 cm[ncm-1][cmuse] = pii ; /* add pii to component */
95 used[pii] = 1 ; /* mark pii as used */
96 }
97 }
98 } /* this component is done */
99
100 cm[ncm-1][0] = cmuse ; /* store size of component */
101
102 if( cmall > cmuse+1 ) /* trim component array */
103 cm[ncm-1] = (int *) realloc( cm[ncm-1], sizeof(int)*(cmuse+1) ) ;
104
105 } while( 1 ) ; /* end of loop over components */
106
107 /* prepare for output */
108
109 free(used) ; /* toss the trash */
110
111 if( ncm == 0 ){ free(cm) ; cm = NULL ; } /* no components!? */
112
113 /* sort components by size */
114
115 for( kk=0 ; kk < ncm ; kk++ ){
116 for( ijk=0,ii=1 ; ii < ncm ; ii++ ){
117 if( cm[ii-1][0] < cm[ii][0] ){
118 tcm = cm[ii-1] ; cm[ii-1] = cm[ii] ; cm[ii] = tcm ; ijk++ ;
119 }
120 }
121 if( !ijk ) break ;
122 }
123
124 *ncom = ncm ; *cmat = cm ; return ;
125 }
126