Doxygen Source Code Documentation
edt_clust.c File Reference
#include "mrilib.h"Go to the source code of this file.
Functions | |
| MCW_cluster_array * | MCW_find_clusters (int nx, int ny, int nz, float dx, float dy, float dz, int ftype, void *fim, float max_dist) |
| void | MCW_cluster_to_vol (int nx, int ny, int nz, int ftype, void *fim, MCW_cluster *clust) |
| void | MCW_erode_clusters (int nx, int ny, int nz, float dx, float dy, float dz, int ftype, void *fim, float max_dist, float pv, int dilate) |
Function Documentation
|
||||||||||||||||||||||||||||
|
Definition at line 200 of file edt_clust.c. References ENTRY, fim, MCW_cluster::i, MCW_cluster::j, MCW_cluster::k, MCW_cluster::mag, MCW_cluster::num_pt, nz, and THREE_TO_IJK. Referenced by EDIT_one_dataset(), main(), RCREND_reload_func_dset(), and REND_reload_func_dset().
00202 {
00203 int icl, ijk ;
00204 int nxy ;
00205 short * sfar ;
00206 float * ffar ;
00207 byte * bfar ;
00208
00209 ENTRY("MCW_cluster_to_vol") ;
00210
00211 if( fim == NULL || clust == NULL ) EXRETURN ;
00212
00213 nxy = nx * ny;
00214
00215 switch( ftype ){
00216 case MRI_short:
00217 sfar = (short *) fim ;
00218 for( icl=0 ; icl < clust->num_pt ; icl++ )
00219 {
00220 ijk = THREE_TO_IJK (clust->i[icl], clust->j[icl], clust->k[icl],
00221 nx, nxy);
00222 sfar[ijk] = clust->mag[icl] ;
00223 }
00224 EXRETURN ;
00225
00226 case MRI_byte:
00227 bfar = (byte *) fim ;
00228 for( icl=0 ; icl < clust->num_pt ; icl++ )
00229 {
00230 ijk = THREE_TO_IJK (clust->i[icl], clust->j[icl], clust->k[icl],
00231 nx, nxy);
00232 bfar[ijk] = clust->mag[icl] ;
00233 }
00234 EXRETURN ;
00235
00236 case MRI_float:
00237 ffar = (float *) fim ;
00238 for( icl=0 ; icl < clust->num_pt ; icl++ )
00239 {
00240 ijk = THREE_TO_IJK (clust->i[icl], clust->j[icl], clust->k[icl],
00241 nx, nxy);
00242 ffar[ijk] = clust->mag[icl] ;
00243 }
00244 EXRETURN ;
00245 }
00246
00247 EXRETURN ; /* should not be reached */
00248 }
|
|
||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 271 of file edt_clust.c. References ENTRY, fim, free, MCW_cluster::i, IJK_TO_THREE, MCW_cluster::j, MCW_cluster::k, KILL_CLUSTER, malloc, MCW_build_mask(), MCW_cluster::num_pt, nz, and THREE_TO_IJK. Referenced by EDIT_one_dataset().
00282 {
00283 MCW_cluster * mask = NULL; /* mask determines nbhd membership */
00284 int minimum; /* minimum number of voxels in nbhd */
00285 int count; /* count of voxels in neighborhood */
00286 int nxy, nxyz; /* numbers of voxels */
00287 int ijk, iv, jv, kv; /* voxel indices */
00288 int ijkm, im, jm, km; /* voxel indices */
00289 int imask, nmask; /* mask indices */
00290 short * sfar; /* pointer to short data */
00291 byte * bfar; /* pointer to byte data */
00292 float * ffar; /* pointer to float data */
00293 float * efim = NULL; /* copy of eroded voxels */
00294
00295 ENTRY("MCW_erode_clusters") ;
00296
00297
00298 /*----- Just in case -----*/
00299 if ( fim == NULL ) EXRETURN;
00300
00301
00302 /*----- Initialize local variables -----*/
00303 nxy = nx * ny; nxyz = nxy * nz;
00304
00305
00306 /*----- Set pointer to input data -----*/
00307 switch (ftype)
00308 {
00309 default: EXRETURN;
00310 case MRI_short: sfar = (short *) fim; break;
00311 case MRI_byte : bfar = (byte *) fim; break;
00312 case MRI_float: ffar = (float *) fim; break;
00313 }
00314
00315
00316 /*----- Initialization for copy of eroded voxels -----*/
00317 efim = (float *) malloc (sizeof(float) * nxyz);
00318 if (efim == NULL)
00319 {
00320 fprintf (stderr, "Unable to allocate memory in MCW_erode_clusters");
00321 EXRETURN;
00322 }
00323 for (ijk = 0; ijk < nxyz; ijk++)
00324 efim[ijk] = 0.0;
00325
00326
00327 /*--- Make a cluster that is a mask of points closer than max_dist ---*/
00328 mask = MCW_build_mask (nx, ny, nz, dx, dy, dz, max_dist);
00329 if (mask == NULL)
00330 {
00331 fprintf (stderr, "Unable to build mask in MCW_erode_clusters");
00332 EXRETURN;
00333 }
00334
00335
00336 /*----- Calculate minimum number of voxels in nbhd. for non-erosion -----*/
00337 nmask = mask->num_pt ;
00338 minimum = floor(pv*nmask + 0.99);
00339 if (minimum <= 0) EXRETURN; /*----- Nothing will be eroded -----*/
00340
00341
00342 /*----- Step 1: Identify voxels to be eroded -----*/
00343 switch (ftype)
00344 {
00345 case MRI_short:
00346 for (ijk = 0; ijk < nxyz; ijk++)
00347 {
00348 if (sfar[ijk] == 0) continue;
00349 IJK_TO_THREE (ijk, iv, jv, kv, nx, nxy);
00350
00351 /*----- Count number of active voxels in the neighborhood -----*/
00352 count = 0;
00353 for (imask = 0; imask < nmask; imask++)
00354 {
00355 im = iv + mask->i[imask];
00356 jm = jv + mask->j[imask];
00357 km = kv + mask->k[imask];
00358 if ( im < 0 || jm < 0 || km < 0 ||
00359 im >= nx || jm >= ny || km >= nz ) continue;
00360 ijkm = THREE_TO_IJK (im, jm, km, nx, nxy);
00361 if (sfar[ijkm] != 0) count++;
00362 }
00363
00364 /*----- Record voxel to be eroded -----*/
00365 if (count < minimum) efim[ijk] = (float) sfar[ijk];
00366 }
00367 break;
00368
00369 case MRI_byte:
00370 for (ijk = 0; ijk < nxyz; ijk++)
00371 {
00372 if (bfar[ijk] == 0) continue;
00373 IJK_TO_THREE (ijk, iv, jv, kv, nx, nxy);
00374
00375 /*----- Count number of active voxels in the neighborhood -----*/
00376 count = 0;
00377 for (imask = 0; imask < nmask; imask++)
00378 {
00379 im = iv + mask->i[imask];
00380 jm = jv + mask->j[imask];
00381 km = kv + mask->k[imask];
00382 if ( im < 0 || jm < 0 || km < 0 ||
00383 im >= nx || jm >= ny || km >= nz ) continue;
00384 ijkm = THREE_TO_IJK (im, jm, km, nx, nxy);
00385 if (bfar[ijkm] != 0) count++;
00386 }
00387
00388 /*----- Record voxel to be eroded -----*/
00389 if (count < minimum) efim[ijk] = (float) bfar[ijk];
00390 }
00391 break;
00392
00393 case MRI_float:
00394 for (ijk = 0; ijk < nxyz; ijk++)
00395 {
00396 if (ffar[ijk] == 0.0) continue;
00397 IJK_TO_THREE (ijk, iv, jv, kv, nx, nxy);
00398
00399 /*----- Count number of active voxels in the neighborhood -----*/
00400 count = 0;
00401 for (imask = 0; imask < nmask; imask++)
00402 {
00403 im = iv + mask->i[imask];
00404 jm = jv + mask->j[imask];
00405 km = kv + mask->k[imask];
00406 if ( im < 0 || jm < 0 || km < 0 ||
00407 im >= nx || jm >= ny || km >= nz ) continue;
00408 ijkm = THREE_TO_IJK (im, jm, km, nx, nxy);
00409 if (ffar[ijkm] != 0.0) count++;
00410 }
00411
00412 /*----- Record voxel to be eroded -----*/
00413 if (count < minimum) efim[ijk] = ffar[ijk];
00414 }
00415 break;
00416
00417 }
00418
00419
00420 /*----- Step 2: Erode voxels -----*/
00421 switch (ftype)
00422 {
00423 case MRI_short:
00424 for (ijk = 0; ijk < nxyz; ijk++)
00425 if (efim[ijk] != 0.0) sfar[ijk] = 0;
00426 break;
00427
00428 case MRI_byte:
00429 for (ijk = 0; ijk < nxyz; ijk++)
00430 if (efim[ijk] != 0.0) bfar[ijk] = 0;
00431 break;
00432
00433 case MRI_float:
00434 for (ijk = 0; ijk < nxyz; ijk++)
00435 if (efim[ijk] != 0.0) ffar[ijk] = 0.0;
00436 break;
00437 }
00438
00439
00440 /*----- Proceed with dilation phase? -----*/
00441 if (dilate)
00442 {
00443
00444
00445 /*----- Step 3: Identify voxels to be restored -----*/
00446 switch (ftype)
00447 {
00448 case MRI_short:
00449 for (ijk = 0; ijk < nxyz; ijk++)
00450 {
00451 if (efim[ijk] == 0.0) continue;
00452 IJK_TO_THREE (ijk, iv, jv, kv, nx, nxy);
00453
00454 /*---- Determine if any active voxels in the neighborhood ----*/
00455 for (imask = 0; imask < nmask; imask++)
00456 {
00457 im = iv + mask->i[imask];
00458 jm = jv + mask->j[imask];
00459 km = kv + mask->k[imask];
00460 if ( im < 0 || jm < 0 || km < 0 ||
00461 im >= nx || jm >= ny || km >= nz ) continue;
00462 ijkm = THREE_TO_IJK (im, jm, km, nx, nxy);
00463 if (sfar[ijkm] != 0) break;
00464 }
00465
00466 /*----- Reset voxel not to be restored -----*/
00467 if (imask == nmask) efim[ijk] = 0.0;
00468 }
00469 break;
00470
00471 case MRI_byte:
00472 for (ijk = 0; ijk < nxyz; ijk++)
00473 {
00474 if (efim[ijk] == 0.0) continue;
00475 IJK_TO_THREE (ijk, iv, jv, kv, nx, nxy);
00476
00477 /*---- Determine if any active voxels in the neighborhood ----*/
00478 for (imask = 0; imask < nmask; imask++)
00479 {
00480 im = iv + mask->i[imask];
00481 jm = jv + mask->j[imask];
00482 km = kv + mask->k[imask];
00483 if ( im < 0 || jm < 0 || km < 0 ||
00484 im >= nx || jm >= ny || km >= nz ) continue;
00485 ijkm = THREE_TO_IJK (im, jm, km, nx, nxy);
00486 if (bfar[ijkm] != 0) break;
00487 }
00488
00489 /*----- Reset voxel not to be restored -----*/
00490 if (imask == nmask) efim[ijk] = 0.0;
00491 }
00492 break;
00493
00494 case MRI_float:
00495 for (ijk = 0; ijk < nxyz; ijk++)
00496 {
00497 if (efim[ijk] == 0.0) continue;
00498 IJK_TO_THREE (ijk, iv, jv, kv, nx, nxy);
00499
00500 /*---- Determine if any active voxels in the neighborhood ----*/
00501 for (imask = 0; imask < nmask; imask++)
00502 {
00503 im = iv + mask->i[imask];
00504 jm = jv + mask->j[imask];
00505 km = kv + mask->k[imask];
00506 if ( im < 0 || jm < 0 || km < 0 ||
00507 im >= nx || jm >= ny || km >= nz ) continue;
00508 ijkm = THREE_TO_IJK (im, jm, km, nx, nxy);
00509 if (ffar[ijkm] != 0.0) break;
00510 }
00511
00512 /*----- Reset voxel not to be restored -----*/
00513 if (imask == nmask) efim[ijk] = 0.0;
00514 }
00515 break;
00516 }
00517
00518
00519 /*----- Step 4: Restore voxels -----*/
00520 switch (ftype)
00521 {
00522 case MRI_short:
00523 for (ijk = 0; ijk < nxyz; ijk++)
00524 if (efim[ijk] != 0.0) sfar[ijk] = (short) efim[ijk];
00525 break;
00526
00527 case MRI_byte:
00528 for (ijk = 0; ijk < nxyz; ijk++)
00529 if (efim[ijk] != 0.0) bfar[ijk] = (byte) efim[ijk];
00530 break;
00531
00532 case MRI_float:
00533 for (ijk = 0; ijk < nxyz; ijk++)
00534 if (efim[ijk] != 0.0) ffar[ijk] = efim[ijk];
00535 break;
00536 }
00537
00538 } /* if (dilate) */
00539
00540
00541 /*----- Release memory -----*/
00542 KILL_CLUSTER(mask) ;
00543 free (efim); efim = NULL;
00544 EXRETURN ;
00545 }
|
|
||||||||||||||||||||||||||||||||||||||||
|
Definition at line 23 of file edt_clust.c. References ADDTO_CLARR, ADDTO_CLUSTER, DESTROY_CLARR, ENTRY, fim, MCW_cluster::i, IJK_TO_THREE, INIT_CLARR, INIT_CLUSTER, MCW_cluster::j, MCW_cluster::k, KILL_CLUSTER, MCW_build_mask(), MCW_cluster_array::num_clu, MCW_cluster::num_pt, nz, RETURN, and THREE_TO_IJK. Referenced by EDIT_one_dataset(), identify_clusters(), main(), NIH_find_clusters(), RCREND_reload_func_dset(), REND_reload_func_dset(), and ROIPLOT_main().
00028 {
00029 MCW_cluster_array * clust_arr ;
00030 MCW_cluster * clust , * mask ;
00031 int ii,jj,kk , nxy,nxyz , ijk , ijk_last , mnum ;
00032 int icl , jma , ijkcl , ijkma , did_one ;
00033 float fimv ;
00034 short * sfar ;
00035 float * ffar ;
00036 byte * bfar ;
00037 short ic, jc, kc;
00038 short im, jm, km;
00039
00040 ENTRY("MCW_find_clusters") ;
00041
00042 if( fim == NULL ) RETURN(NULL) ;
00043
00044 switch( ftype ){
00045 default: RETURN(NULL) ;
00046 case MRI_short: sfar = (short *) fim ; break ;
00047 case MRI_byte : bfar = (byte *) fim ; break ;
00048 case MRI_float: ffar = (float *) fim ; break ;
00049 }
00050
00051 /*--- make a cluster that is a mask of points closer than max_dist ---*/
00052
00053 mask = MCW_build_mask (nx, ny, nz, dx, dy, dz, max_dist);
00054 if (mask == NULL)
00055 {
00056 fprintf (stderr, "Unable to build mask in MCW_find_clusters");
00057 RETURN(NULL) ;
00058 }
00059
00060 nxy = nx*ny ; nxyz = nxy * nz ;
00061
00062 mnum = mask->num_pt ;
00063
00064
00065 /*--- scan through array, find nonzero point, build a cluster, ... ---*/
00066
00067 INIT_CLARR(clust_arr) ;
00068
00069 ijk_last = 0 ;
00070 do {
00071 switch( ftype ){
00072 case MRI_short:
00073 for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( sfar[ijk] != 0 ) break ;
00074 if( ijk < nxyz ){
00075 fimv = sfar[ijk] ; sfar[ijk] = 0 ; /* save found point */
00076 }
00077 break ;
00078
00079 case MRI_byte:
00080 for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( bfar[ijk] != 0 ) break ;
00081 if( ijk < nxyz ){
00082 fimv = bfar[ijk] ; bfar[ijk] = 0 ; /* save found point */
00083 }
00084 break ;
00085
00086 case MRI_float:
00087 for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( ffar[ijk] != 0.0 ) break ;
00088 if( ijk < nxyz ){
00089 fimv = ffar[ijk] ; ffar[ijk] = 0.0 ; /* save found point */
00090 }
00091 break ;
00092 }
00093 if( ijk == nxyz ) break ; /* didn't find any! */
00094
00095 #ifdef CLUST_DEBUG
00096 printf(" starting cluster at ijk=%d\n",ijk) ;
00097 #endif
00098
00099 ijk_last = ijk+1 ; /* start here next time */
00100
00101 INIT_CLUSTER(clust) ; /* make a new cluster */
00102 IJK_TO_THREE(ijk,ic,jc,kc,nx,nxy) ;
00103 ADDTO_CLUSTER( clust , ic, jc, kc, fimv ) ; /* start it off */
00104
00105 /*--
00106 for each point in cluster:
00107 check points offset by the mask for nonzero entries in fim
00108 enter those into cluster
00109 continue until end of cluster is reached
00110 (note that cluster is expanding as we progress)
00111 --*/
00112
00113 switch( ftype ){
00114 case MRI_short:
00115 for( icl=0 ; icl < clust->num_pt ; icl++ ){
00116 ic = clust->i[icl];
00117 jc = clust->j[icl];
00118 kc = clust->k[icl];
00119
00120 for( jma=0 ; jma < mnum ; jma++ ){
00121 im = ic + mask->i[jma];
00122 jm = jc + mask->j[jma];
00123 km = kc + mask->k[jma];
00124 if( im < 0 || im >= nx ||
00125 jm < 0 || jm >= ny || km < 0 || km >= nz ) continue ;
00126
00127 ijkma = THREE_TO_IJK (im, jm, km, nx, nxy);
00128 if( ijkma < ijk_last || ijkma >= nxyz || sfar[ijkma] == 0 ) continue ;
00129
00130 ADDTO_CLUSTER( clust , im, jm, km, sfar[ijkma] ) ;
00131 sfar[ijkma] = 0 ;
00132 }
00133 }
00134 break ;
00135
00136 case MRI_byte:
00137 for( icl=0 ; icl < clust->num_pt ; icl++ ){
00138 ic = clust->i[icl];
00139 jc = clust->j[icl];
00140 kc = clust->k[icl];
00141
00142 for( jma=0 ; jma < mnum ; jma++ ){
00143 im = ic + mask->i[jma];
00144 jm = jc + mask->j[jma];
00145 km = kc + mask->k[jma];
00146 if( im < 0 || im >= nx ||
00147 jm < 0 || jm >= ny || km < 0 || km >= nz ) continue ;
00148
00149 ijkma = THREE_TO_IJK (im, jm, km, nx, nxy);
00150 if( ijkma < ijk_last || ijkma >= nxyz || bfar[ijkma] == 0 ) continue ;
00151
00152 ADDTO_CLUSTER( clust , im, jm, km, bfar[ijkma] ) ;
00153 bfar[ijkma] = 0 ;
00154 }
00155 }
00156 break ;
00157
00158 case MRI_float:
00159 for( icl=0 ; icl < clust->num_pt ; icl++ ){
00160 ic = clust->i[icl];
00161 jc = clust->j[icl];
00162 kc = clust->k[icl];
00163
00164 for( jma=0 ; jma < mnum ; jma++ ){
00165 im = ic + mask->i[jma];
00166 jm = jc + mask->j[jma];
00167 km = kc + mask->k[jma];
00168 if( im < 0 || im >= nx ||
00169 jm < 0 || jm >= ny || km < 0 || km >= nz ) continue ;
00170
00171 ijkma = THREE_TO_IJK (im, jm, km, nx, nxy);
00172 if( ijkma < ijk_last || ijkma >= nxyz || ffar[ijkma] == 0.0 ) continue ;
00173
00174 ADDTO_CLUSTER( clust , im, jm, km, ffar[ijkma] ) ;
00175 ffar[ijkma] = 0.0 ;
00176 }
00177 }
00178 break ;
00179 }
00180
00181 ADDTO_CLARR(clust_arr,clust) ;
00182 } while( 1 ) ;
00183
00184 KILL_CLUSTER(mask) ;
00185
00186 if( clust_arr->num_clu <= 0 ){ DESTROY_CLARR(clust_arr) ; }
00187
00188 RETURN(clust_arr) ;
00189 }
|
