00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 MCW_cluster_array * MCW_find_clusters(
00024 int nx, int ny, int nz,
00025 float dx, float dy, float dz,
00026 int ftype , void * fim ,
00027 float max_dist )
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
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
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 ;
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 ;
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 ;
00090 }
00091 break ;
00092 }
00093 if( ijk == nxyz ) break ;
00094
00095 #ifdef CLUST_DEBUG
00096 printf(" starting cluster at ijk=%d\n",ijk) ;
00097 #endif
00098
00099 ijk_last = ijk+1 ;
00100
00101 INIT_CLUSTER(clust) ;
00102 IJK_TO_THREE(ijk,ic,jc,kc,nx,nxy) ;
00103 ADDTO_CLUSTER( clust , ic, jc, kc, fimv ) ;
00104
00105
00106
00107
00108
00109
00110
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 }
00190
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200 void MCW_cluster_to_vol( int nx , int ny , int nz ,
00201 int ftype , void * fim , MCW_cluster * clust )
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 ;
00248 }
00249
00250
00251
00252
00253
00254
00255
00256
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268
00269
00270 void MCW_erode_clusters
00271 (
00272 int nx, int ny, int nz,
00273 float dx, float dy, float dz,
00274 int ftype,
00275 void * fim,
00276 float max_dist,
00277 float pv,
00278
00279 int dilate
00280 )
00281
00282 {
00283 MCW_cluster * mask = NULL;
00284 int minimum;
00285 int count;
00286 int nxy, nxyz;
00287 int ijk, iv, jv, kv;
00288 int ijkm, im, jm, km;
00289 int imask, nmask;
00290 short * sfar;
00291 byte * bfar;
00292 float * ffar;
00293 float * efim = NULL;
00294
00295 ENTRY("MCW_erode_clusters") ;
00296
00297
00298
00299 if ( fim == NULL ) EXRETURN;
00300
00301
00302
00303 nxy = nx * ny; nxyz = nxy * nz;
00304
00305
00306
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
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
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
00337 nmask = mask->num_pt ;
00338 minimum = floor(pv*nmask + 0.99);
00339 if (minimum <= 0) EXRETURN;
00340
00341
00342
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
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
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
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
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
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
00413 if (count < minimum) efim[ijk] = ffar[ijk];
00414 }
00415 break;
00416
00417 }
00418
00419
00420
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
00441 if (dilate)
00442 {
00443
00444
00445
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
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
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
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
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
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
00513 if (imask == nmask) efim[ijk] = 0.0;
00514 }
00515 break;
00516 }
00517
00518
00519
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 }
00539
00540
00541
00542 KILL_CLUSTER(mask) ;
00543 free (efim); efim = NULL;
00544 EXRETURN ;
00545 }