1 #include "mrilib.h"
2
3 /*---------------------------------------------------------------------*/
4 static char *report = NULL ;
mri_clusterize_report(void)5 char * mri_clusterize_report(void){ return report; }
6 /*---------------------------------------------------------------------*/
7
8 #if 0
9 MCW_cluster_array *find_clusters_NN1( MRI_IMAGE *cim ) ;
10 MCW_cluster_array *find_clusters_NN2( MRI_IMAGE *cim ) ;
11 MCW_cluster_array *find_clusters_NN3( MRI_IMAGE *cim ) ;
12 #endif
13
14 /*---------------------------------------------------------------------*/
15 static MCW_cluster_array *clarout=NULL ;
16
mri_clusterize_array(int clear)17 MCW_cluster_array * mri_clusterize_array(int clear)
18 {
19 MCW_cluster_array *cc = clarout ;
20 if( clear ) clarout = NULL ;
21 return cc ;
22 }
23
24 /*---------------------------------------------------------------------*/
25 /*! Cluster-edit volume bim, possibly thresholding with tim, and
26 produce a new output image. [05 Sep 2006]
27
28 Modified [12 Jul 2021] to allow 1-voxel 'clusters' (for DRG).
29 -----------------------------------------------------------------------*/
30
mri_clusterize(float rmm,float vmul,MRI_IMAGE * bim,float thb,float tht,MRI_IMAGE * tim,int posonly,byte * mask)31 MRI_IMAGE * mri_clusterize( float rmm , float vmul , MRI_IMAGE *bim ,
32 float thb , float tht , MRI_IMAGE *tim ,
33 int posonly , byte *mask )
34 {
35 float dx,dy,dz , dbot ;
36 int nx,ny,nz , ptmin,iclu , nkeep,nkill,ncgood , nbot,ntop , ii ;
37 MRI_IMAGE *cim ; void *car ;
38 MCW_cluster *cl , *dl ; MCW_cluster_array *clar ;
39 int nnlev = 0 ;
40 static char *cclev[3] = { "[faces touch]" ,
41 "[edges touch]" , "[corners touch]" } ;
42
43 ENTRY("mri_clusterize") ;
44
45 if( report != NULL ){ free((void *)report); report = NULL; }
46 if( clarout != NULL ){ DESTROY_CLARR(clarout); }
47
48 if( bim == NULL || mri_data_pointer(bim) == NULL ) RETURN(NULL) ;
49
50 nx = bim->nx; ny = bim->ny; nz = bim->nz;
51 dx = bim->dx; dy = bim->dy; dz = bim->dz;
52 if( dx <= 0.0f ) dx = 1.0f ;
53 if( dy <= 0.0f ) dy = 1.0f ;
54 if( dz <= 0.0f ) dz = 1.0f ;
55 dbot = MIN(dx,dy) ; dbot = MIN(dbot,dz) ;
56
57 if( rmm < dbot ){
58 int irr = (int)rintf(rmm) ;
59 dx = dy = dz = 1.0f;
60 switch( irr ){
61 default: rmm = 1.01f ; nnlev = 1 ;break ; /* NN1 */
62 case -2: rmm = 1.44f ; nnlev = 2 ;break ; /* NN2 */
63 case -3: rmm = 1.75f ; nnlev = 3 ;break ; /* NN3 */
64 }
65 }
66
67 /* create copy of input image (this will be edited below) */
68
69 cim = mri_copy(bim) ; car = mri_data_pointer(cim) ;
70 if( car == NULL ){ mri_free(cim) ; RETURN(NULL) ; }
71
72 /* threshold it, if so ordered (note that tim==bim is legal) */
73
74 if( tht > thb && tim != NULL )
75 mri_threshold( thb , tht , tim , cim ) ;
76
77 if( posonly )
78 mri_threshold( -1.e9 , 0.0 , cim , cim ) ;
79
80 mri_maskify( cim , mask ) ; /* Jul 2010: mask it? */
81
82 /* smallest cluster to keep */
83
84 ptmin = (int)( vmul/(dx*dy*dz) + 0.99f ) ;
85 if( ptmin < 1 ) ptmin = 1 ;
86
87 /* find all clusters */
88
89 clar = MCW_find_clusters( nx,ny,nz , dx,dy,dz , cim->kind,car , rmm ) ;
90
91 /* put all big clusters back into the image */
92
93 if( clar != NULL ){
94 ncgood = nkeep = nkill = 0 ; nbot = 999999 ; ntop = 0 ;
95 for( iclu=0 ; iclu < clar->num_clu ; iclu++ ){
96 cl = clar->clar[iclu] ;
97 if( cl->num_pt >= ptmin ){ /* put back into array, get some stats */
98
99 MCW_cluster_to_vol( nx,ny,nz , cim->kind,car , cl ) ;
100
101 nkeep += cl->num_pt ; ncgood++ ;
102 nbot = MIN(nbot,cl->num_pt); ntop = MAX(ntop,cl->num_pt);
103
104 if( clarout == NULL ) INIT_CLARR(clarout) ;
105 COPY_CLUSTER(dl,cl) ; ADDTO_CLARR(clarout,dl) ;
106 } else {
107 nkill += cl->num_pt ;
108 }
109 }
110 DESTROY_CLARR(clar) ;
111 report = THD_zzprintf( report ,
112 " Voxels survived clustering =%6d\n"
113 " Voxels edited out =%6d\n" ,
114 nkeep , nkill ) ;
115 #if 0
116 if( ntop >= nbot )
117 report = THD_zzprintf( report ,
118 " Min cluster size (voxels) =%6d\n"
119 " Max cluster size =%6d\n"
120 " Number of clusters kept =%6d\n" ,
121 nbot , ntop , ncgood ) ;
122 #endif
123 if( nnlev > 0 )
124 report = THD_zzprintf( report ,
125 " NN clustering level =%6d %s\n" ,
126 nnlev , cclev[nnlev-1] ) ;
127 }
128
129 RETURN(cim) ;
130 }
131
132 /*---------------------------------------------------------------------*/
133 /*! Bi-Cluster-edit volume bim, possibly thresholding with tim, and
134 produce a new output image. [05 Sep 2006]
135 -----------------------------------------------------------------------*/
136
mri_bi_clusterize(float rmm,float vmul,MRI_IMAGE * bim,float thb,float tht,MRI_IMAGE * tim,byte * mask)137 MRI_IMAGE * mri_bi_clusterize( float rmm , float vmul , MRI_IMAGE *bim ,
138 float thb , float tht , MRI_IMAGE *tim ,
139 byte *mask )
140 {
141 float dx,dy,dz , dbot ;
142 int nx,ny,nz , ptmin,iclu , nkeep=0,nkill=0,ncgood=0 ;
143 int nbot=9999999,ntop=0 , ii , pclust=0,nclust=0 ;
144 MRI_IMAGE *cim , *dim ; void *car , *dar ;
145 MCW_cluster *cl , *dl ; MCW_cluster_array *clar ;
146 int nnlev = 0 ;
147 static char *cclev[3] = { "[faces touch]" ,
148 "[edges touch]" , "[corners touch]" } ;
149
150 ENTRY("mri_bi_clusterize") ;
151
152 if( report != NULL ){ free((void *)report); report = NULL; }
153 if( clarout != NULL ){ DESTROY_CLARR(clarout); }
154
155 if( bim == NULL || mri_data_pointer(bim) == NULL ) RETURN(NULL) ;
156
157 if( thb > 0.0f || tht < 0.0f || tim == NULL ) RETURN(NULL) ;
158
159 nx = bim->nx; ny = bim->ny; nz = bim->nz;
160 dx = bim->dx; dy = bim->dy; dz = bim->dz;
161 if( dx <= 0.0f ) dx = 1.0f ;
162 if( dy <= 0.0f ) dy = 1.0f ;
163 if( dz <= 0.0f ) dz = 1.0f ;
164 dbot = MIN(dx,dy) ; dbot = MIN(dbot,dz) ;
165
166 if( rmm < dbot ){
167 int irr = (int)rintf(rmm) ;
168 dx = dy = dz = 1.0f;
169 switch( irr ){
170 default: rmm = 1.01f ; nnlev = 1 ;break ; /* NN1 */
171 case -2: rmm = 1.44f ; nnlev = 2 ;break ; /* NN2 */
172 case -3: rmm = 1.75f ; nnlev = 3 ;break ; /* NN3 */
173 }
174 }
175 ptmin = (int)( vmul/(dx*dy*dz) + 0.99f ) ; /* smallest cluster to keep */
176 if( ptmin < 1 ) ptmin = 1 ;
177
178 /* 0-filled copy of input == will be output image */
179
180 dim = mri_new_conforming(bim,bim->kind) ; dar = mri_data_pointer(dim) ;
181
182 /* create copy of input image to be edited below */
183
184 cim = mri_copy(bim) ; car = mri_data_pointer(cim) ;
185 mri_maskify( cim , mask ) ; /* mask it? */
186
187 /* threshold it on the positive side */
188
189 mri_threshold( -1.e22 , tht , tim , cim ) ;
190
191 /* find all positive clusters */
192
193 clar = MCW_find_clusters( nx,ny,nz , dx,dy,dz , cim->kind,car , rmm ) ;
194 mri_free(cim) ;
195
196 /* put all big clusters into the output image */
197
198 if( clar != NULL ){
199 for( iclu=0 ; iclu < clar->num_clu ; iclu++ ){
200 cl = clar->clar[iclu] ;
201 if( cl->num_pt >= ptmin ){ /* put back into array, get some stats */
202
203 MCW_cluster_to_vol( nx,ny,nz , dim->kind,dar , cl ) ;
204
205 nkeep += cl->num_pt ; ncgood++ ; pclust++ ;
206 nbot = MIN(nbot,cl->num_pt); ntop = MAX(ntop,cl->num_pt);
207
208 if( clarout == NULL ) INIT_CLARR(clarout) ;
209 COPY_CLUSTER(dl,cl) ; ADDTO_CLARR(clarout,dl) ;
210 } else {
211 nkill += cl->num_pt ;
212 }
213 }
214 DESTROY_CLARR(clar) ;
215 }
216
217 /* repeat for negative clusters */
218
219 cim = mri_copy(bim) ; car = mri_data_pointer(cim) ;
220 mri_maskify( cim , mask ) ;
221 mri_threshold( thb , 1.e+22 , tim , cim ) ;
222 clar = MCW_find_clusters( nx,ny,nz , dx,dy,dz , cim->kind,car , rmm ) ;
223 mri_free(cim) ;
224
225 if( clar != NULL ){
226 for( iclu=0 ; iclu < clar->num_clu ; iclu++ ){
227 cl = clar->clar[iclu] ;
228 if( cl->num_pt >= ptmin ){ /* put back into array, get some stats */
229 MCW_cluster_to_vol( nx,ny,nz , dim->kind,dar , cl ) ;
230 nkeep += cl->num_pt ; ncgood++ ; nclust++ ;
231 nbot = MIN(nbot,cl->num_pt); ntop = MAX(ntop,cl->num_pt);
232 if( clarout == NULL ) INIT_CLARR(clarout) ;
233 COPY_CLUSTER(dl,cl) ; ADDTO_CLARR(clarout,dl) ;
234 } else {
235 nkill += cl->num_pt ;
236 }
237 }
238 DESTROY_CLARR(clar) ;
239 }
240
241 if( nkeep > 0 ){
242 report = THD_zzprintf( report ,
243 " Voxels survived clustering =%6d\n"
244 " Voxels edited out =%6d\n"
245 " # Positive clusters =%6d\n"
246 " # Negative clusters =%6d\n" ,
247 nkeep , nkill , pclust,nclust ) ;
248 if( nnlev > 0 )
249 report = THD_zzprintf( report ,
250 " NN clustering level =%6d %s\n" ,
251 nnlev , cclev[nnlev-1] ) ;
252 }
253
254 RETURN(dim) ;
255 }
256
257 /*---------------------------------------------------------------------*/
258
mri_clusterize_detailize(MCW_cluster * cl,int icent_flag)259 mri_cluster_detail mri_clusterize_detailize( MCW_cluster *cl, int icent_flag )
260 {
261 mri_cluster_detail cld ;
262 float xcm,ycm,zcm , xpk,ypk,zpk , vpk,vvv,vsum ;
263 float xmi,ymi,zmi, xqq,yqq,zqq , wbest ;
264 int ii , kk ;
265
266 ENTRY("mri_clusterize_detailize") ;
267
268 memset( &cld , 0 , sizeof(mri_cluster_detail) ) ;
269 if( cl == NULL || cl->num_pt <= 0 ) RETURN(cld) ;
270
271 cld.nvox = cl->num_pt ;
272 cld.volume = cl->num_pt ;
273 xcm = ycm = zcm = 0.0f ; xpk = ypk = zpk = vpk = 0.0f ;
274 xmi = ymi = zmi = 0.0f ;
275
276 for( vsum=ii=0 ; ii < cl->num_pt ; ii++ ){
277 vvv = fabsf(cl->mag[ii]) ; vsum += vvv ;
278 xcm += vvv*cl->i[ii] ; ycm += vvv*cl->j[ii] ; zcm += vvv*cl->k[ii] ;
279 if( vvv > vpk ){
280 xpk = cl->i[ii]; ypk = cl->j[ii]; zpk = cl->k[ii]; vpk = vvv;
281 }
282 }
283 if( vsum > 0.0f ){
284 xcm = cld.xcm = xcm / vsum;
285 ycm = cld.ycm = ycm / vsum;
286 zcm = cld.zcm = zcm / vsum;
287 }
288 cld.xpk = xpk; cld.ypk = ypk; cld.zpk = zpk;
289
290 /* return early if icent not needed - avoid slow process below*/
291 if( !icent_flag ){
292 cld.xmi = xmi ; cld.ymi = ymi ; cld.zmi = zmi ;
293 RETURN(cld);
294 }
295
296 /* find internal center [08 May 2019] */
297 if( vsum > 0.0f ){
298 wbest = 1.e+37f ; xmi=ymi=zmi = 0.0f ;
299 for( kk=0 ; kk < cl->num_pt ; kk++ ){
300 xqq = cl->i[kk] ; yqq = cl->j[kk] ; zqq = cl->k[kk] ;
301 vsum = SQR(xqq-xcm) + SQR(yqq-ycm) + SQR(zqq-zcm) ;
302 if( vsum < wbest ){ xmi = xqq; ymi = yqq; zmi = zqq; wbest=vsum; }
303 }
304 }
305 cld.xmi = xmi ; cld.ymi = ymi ; cld.zmi = zmi ;
306
307 RETURN(cld) ;
308 }
309
310 #if 0
311 /*---------------------------------------------------------------------------*/
312 /** The stuff below is added Jul 2010, for the new AFNI Clusterize methods. **/
313 /*---------------------------------------------------------------------------*/
314
315 /* Put (i,j,k) into the current cluster, maybe */
316
317 #undef CPUT
318 #define CPUT(i,j,k) \
319 do{ ijk = THREE_TO_IJK(i,j,k,nx,nxy) ; \
320 if( car[ijk] != 0.0f ){ \
321 ADDTO_CLUSTER(clust,i,j,k,car[ijk]) ; \
322 car[ijk] = 0.0f ; \
323 } } while(0)
324
325 /*----------------------------------------------------------------------------*/
326
327 MCW_cluster_array *find_clusters_NN1( MRI_IMAGE *cim )
328 {
329 MCW_cluster_array *clust_arr ;
330 MCW_cluster *clust ;
331 int ii,jj,kk, icl , ijk , ijk_last ;
332 int ip,jp,kp , im,jm,km ;
333 int nx,ny,nz,nxy,nxyz ;
334 float *car = MRI_FLOAT_PTR(cim) ;
335
336 nx = cim->nx ; ny = cim->ny ; cim->nz ; nxy = nx*ny ; nxyz = nxy*nz ;
337
338 INIT_CLARR(clust_arr) ;
339
340 ijk_last = 0 ; /* start scanning at the start */
341
342 while(1) {
343 /* find next nonzero point in mmm array */
344
345 for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( car[ijk] != 0.0f ) break ;
346 if( ijk == nxyz ) break ; /* didn't find any! */
347 ijk_last = ijk+1 ; /* start here next time */
348
349 INIT_CLUSTER(clust) ;
350 IJK_TO_THREE(ijk,ii,jj,kk,nx,nxy) ;
351 ADDTO_CLUSTER( clust , ii,jj,kk , car[ijk] ) ;
352 car[ijk] = 0.0f ;
353
354 /* loop over points in cluster, checking their neighbors,
355 growing the cluster if we find any that belong therein */
356
357 for( icl=0 ; icl < clust->num_pt ; icl++ ){
358 ii = clust->i[icl] ; jj = clust->j[icl] ; kk = clust->k[icl] ;
359 im = ii-1 ; jm = jj-1 ; km = kk-1 ;
360 ip = ii+1 ; jp = jj+1 ; kp = kk+1 ;
361
362 if( im >= 0 ) CPUT(im,jj,kk) ;
363 if( ip < nx ) CPUT(ip,jj,kk) ;
364 if( jm >= 0 ) CPUT(ii,jm,kk) ;
365 if( jp < ny ) CPUT(ii,jp,kk) ;
366 if( km >= 0 ) CPUT(ii,jj,km) ;
367 if( kp < nz ) CPUT(ii,jj,kp) ;
368 }
369 ADDTO_CLARR(clust_arr,clust) ;
370 }
371
372 if( clust_arr->num_clu <= 0 ){ DESTROY_CLARR(clust_arr) ; }
373 return clust_arr ;
374 }
375
376 /*----------------------------------------------------------------------------*/
377
378 MCW_cluster_array *find_clusters_NN2( MRI_IMAGE *cim )
379 {
380 MCW_cluster_array *clust_arr ;
381 MCW_cluster *clust ;
382 int ii,jj,kk, icl , ijk , ijk_last ;
383 int ip,jp,kp , im,jm,km ;
384 int nx,ny,nz,nxy,nxyz ;
385 float *car = MRI_FLOAT_PTR(cim) ;
386
387 nx = cim->nx ; ny = cim->ny ; cim->nz ; nxy = nx*ny ; nxyz = nxy*nz ;
388
389 INIT_CLARR(clust_arr) ;
390
391 ijk_last = 0 ; /* start scanning at the start */
392
393 while(1) {
394 /* find next nonzero point in mmm array */
395
396 for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( car[ijk] != 0.0f ) break ;
397 if( ijk == nxyz ) break ; /* didn't find any! */
398 ijk_last = ijk+1 ; /* start here next time */
399
400 INIT_CLUSTER(clust) ;
401 IJK_TO_THREE(ijk,ii,jj,kk,nx,nxy) ;
402 ADDTO_CLUSTER( clust , ii,jj,kk , car[ijk] ) ;
403 car[ijk] = 0.0f ;
404
405 /* loop over points in cluster, checking their neighbors,
406 growing the cluster if we find any that belong therein */
407
408 for( icl=0 ; icl < clust->num_pt ; icl++ ){
409 ii = clust->i[icl] ; jj = clust->j[icl] ; kk = clust->k[icl] ;
410 im = ii-1 ; jm = jj-1 ; km = kk-1 ;
411 ip = ii+1 ; jp = jj+1 ; kp = kk+1 ;
412
413 if( im >= 0 ){ CPUT(im,jj,kk) ;
414 if( jm >= 0 ) CPUT(im,jm,kk) ; /* 2NN */
415 if( jp < nx ) CPUT(im,jp,kk) ; /* 2NN */
416 if( km >= 0 ) CPUT(im,jj,km) ; /* 2NN */
417 if( kp < nz ) CPUT(im,jj,kp) ; /* 2NN */
418 }
419 if( ip < nx ){ CPUT(ip,jj,kk) ;
420 if( jm >= 0 ) CPUT(ip,jm,kk) ; /* 2NN */
421 if( jp < nx ) CPUT(ip,jp,kk) ; /* 2NN */
422 if( km >= 0 ) CPUT(ip,jj,km) ; /* 2NN */
423 if( kp < nz ) CPUT(ip,jj,kp) ; /* 2NN */
424 }
425 if( jm >= 0 ){ CPUT(ii,jm,kk) ;
426 if( km >= 0 ) CPUT(ii,jm,km) ; /* 2NN */
427 if( kp < nz ) CPUT(ii,jm,kp) ; /* 2NN */
428 }
429 if( jp < ny ){ CPUT(ii,jp,kk) ;
430 if( km >= 0 ) CPUT(ii,jp,km) ; /* 2NN */
431 if( kp < nz ) CPUT(ii,jp,kp) ; /* 2NN */
432 }
433 if( km >= 0 ) CPUT(ii,jj,km) ;
434 if( kp < nz ) CPUT(ii,jj,kp) ;
435 }
436 ADDTO_CLARR(clust_arr,clust) ;
437 }
438
439 if( clust_arr->num_clu <= 0 ){ DESTROY_CLARR(clust_arr) ; }
440 return clust_arr ;
441 }
442
443 /*----------------------------------------------------------------------------*/
444
445 MCW_cluster_array *find_clusters_NN3( MRI_IMAGE *cim )
446 {
447 MCW_cluster_array *clust_arr ;
448 MCW_cluster *clust ;
449 int ii,jj,kk, icl , ijk , ijk_last ;
450 int ip,jp,kp , im,jm,km ;
451 int nx,ny,nz,nxy,nxyz ;
452 float *car = MRI_FLOAT_PTR(cim) ;
453
454 nx = cim->nx ; ny = cim->ny ; cim->nz ; nxy = nx*ny ; nxyz = nxy*nz ;
455
456 INIT_CLARR(clust_arr) ;
457
458 ijk_last = 0 ; /* start scanning at the start */
459
460 while(1) {
461 /* find next nonzero point in mmm array */
462
463 for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( car[ijk] != 0.0f ) break ;
464 if( ijk == nxyz ) break ; /* didn't find any! */
465 ijk_last = ijk+1 ; /* start here next time */
466
467 INIT_CLUSTER(clust) ;
468 IJK_TO_THREE(ijk,ii,jj,kk,nx,nxy) ;
469 ADDTO_CLUSTER( clust , ii,jj,kk , car[ijk] ) ;
470 car[ijk] = 0.0f ;
471
472 /* loop over points in cluster, checking their neighbors,
473 growing the cluster if we find any that belong therein */
474
475 for( icl=0 ; icl < clust->num_pt ; icl++ ){
476 ii = clust->i[icl] ; jj = clust->j[icl] ; kk = clust->k[icl] ;
477 im = ii-1 ; jm = jj-1 ; km = kk-1 ;
478 ip = ii+1 ; jp = jj+1 ; kp = kk+1 ;
479
480 if( im >= 0 ){ CPUT(im,jj,kk) ;
481 if( jm >= 0 ) CPUT(im,jm,kk) ; /* 2NN */
482 if( jp < nx ) CPUT(im,jp,kk) ; /* 2NN */
483 if( km >= 0 ) CPUT(im,jj,km) ; /* 2NN */
484 if( kp < nz ) CPUT(im,jj,kp) ; /* 2NN */
485 if( jm >= 0 && km >= 0 ) CPUT(im,jm,km) ; /* 3NN */
486 if( jm >= 0 && kp < nz ) CPUT(im,jm,kp) ; /* 3NN */
487 if( jp < ny && km >= 0 ) CPUT(im,jp,km) ; /* 3NN */
488 if( jp < ny && kp < nz ) CPUT(im,jp,kp) ; /* 3NN */
489 }
490 if( ip < nx ){ CPUT(ip,jj,kk) ;
491 if( jm >= 0 ) CPUT(ip,jm,kk) ; /* 2NN */
492 if( jp < nx ) CPUT(ip,jp,kk) ; /* 2NN */
493 if( km >= 0 ) CPUT(ip,jj,km) ; /* 2NN */
494 if( kp < nz ) CPUT(ip,jj,kp) ; /* 2NN */
495 if( jm >= 0 && km >= 0 ) CPUT(ip,jm,km) ; /* 3NN */
496 if( jm >= 0 && kp < nz ) CPUT(ip,jm,kp) ; /* 3NN */
497 if( jp < ny && km >= 0 ) CPUT(ip,jp,km) ; /* 3NN */
498 if( jp < ny && kp < nz ) CPUT(ip,jp,kp) ; /* 3NN */
499 }
500 if( jm >= 0 ){ CPUT(ii,jm,kk) ;
501 if( km >= 0 ) CPUT(ii,jm,km) ; /* 2NN */
502 if( kp < nz ) CPUT(ii,jm,kp) ; /* 2NN */
503 }
504 if( jp < ny ){ CPUT(ii,jp,kk) ;
505 if( km >= 0 ) CPUT(ii,jp,km) ; /* 2NN */
506 if( kp < nz ) CPUT(ii,jp,kp) ; /* 2NN */
507 }
508 if( km >= 0 ) CPUT(ii,jj,km) ;
509 if( kp < nz ) CPUT(ii,jj,kp) ;
510
511 }
512 ADDTO_CLARR(clust_arr,clust) ;
513 }
514
515 if( clust_arr->num_clu <= 0 ){ DESTROY_CLARR(clust_arr) ; }
516 return clust_arr ;
517 }
518 #endif
519