00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009 void find_base_value( int nxyz , short * sfim , int * base , int * peak ) ;
00010 void remove_isolated_stuff( THD_3dim_dataset * qset ) ;
00011 void xyz_to_ijk( THD_3dim_dataset * ds , float x , float y , float z ,
00012 int * i , int * j , int * k ) ;
00013 void DRAW_2dfiller( int nx , int ny , int ix , int jy , byte * ar ) ;
00014
00015
00016
00017 #define DD(i,j,k) dar[(i)+(j)*nx+(k)*nxy]
00018 #define QQ(i,j,k) qar[(i)+(j)*nx+(k)*nxy]
00019 #define BB(i,j) bar[(i)+(j)*nx]
00020
00021 static int verbose = 0 ;
00022
00023 int main( int argc , char * argv[] )
00024 {
00025 THD_3dim_dataset * dset , * qset ;
00026 short * dar , * qar ;
00027
00028 char prefix[128] = "strip" ;
00029 int base = 0 , hyval = 0 , hyper = 0 ;
00030 double metric_thresh = 1.05 , area_thresh = 9999.0 ;
00031
00032 int iarg , nx,ny,nz,nxy,nxyz , dx,dy,dz , ii,jj,kk, ic,jc,kc , ktop , cc,pp ;
00033 byte * bar ;
00034 double qsum , sum , rr , metric,area ;
00035 int nsum , didit ;
00036
00037
00038
00039 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00040 printf(
00041 "Usage: 3dstrip [options] dataset\n"
00042 "Attempts to strip the scalp tissue from the input anatomical dataset.\n"
00043 "This dataset must be in AC-PC or Talairach coordinates, must be stored\n"
00044 "as shorts, and must have cubical voxels.\n"
00045 "\n"
00046 "Options:\n"
00047 " -base bb = Set all voxels at or below value 'bb' to zero.\n"
00048 " [default = computed from the data histogram]\n"
00049 " -metric mm = Set the metric threshold to 'mm'; the metric is\n"
00050 " used to determine the annularity of the removed\n"
00051 " tissue in each axial slice.\n"
00052 " [default = 1.05 (unitless); must be > 1.0]\n"
00053 " -area aa = Set the area threshold to 'aa' square millimeters;\n"
00054 " if the amount of removed tissue would be larger \n"
00055 " than 'aa', or the metric is larger than 'mm',\n"
00056 " then nothing will be removed in or below the slice.\n"
00057 " [default = 9999]\n"
00058 " -hyper = Flags the elimination of hyperintensity voxels after\n"
00059 " the stripping procedure.\n"
00060 " [default = do not remove these voxels]\n"
00061 " -hyclip = Flags the setting of hyperintensity voxels back to\n"
00062 " the hyperintensity threshold (rather than to zero).\n"
00063 " -hyval hh = Set the hyperintensity threshold to 'hh'; all voxels\n"
00064 " with intensity 'hh' or above will be set to zero\n"
00065 " (for -hyper) or set to this value (for -hyclip).\n"
00066 " [default = computed from the data histogram]\n"
00067 " -prefix pp = Use 'pp' for the prefix of the output dataset.\n"
00068 " [default = 'strip']\n"
00069 " -verbose = Print out progress reports, including the area\n"
00070 " and metric for each axial slice\n"
00071 "\n"
00072 "The metric is average(r**2)/average(r)**2 for all candidate removal\n"
00073 "voxels, where r is measured from the center of the slice. This ratio\n"
00074 "is 1 for a circle, and will be nearly 1 for a thin annular region.\n"
00075 "The default threshold of 1.05 is chosen from experience. The default\n"
00076 "area threshold is also chosen from experience. The goal of the thresholds\n"
00077 "is to avoid removing brain tissue.\n"
00078 ) ;
00079 exit(0) ;
00080 }
00081
00082
00083
00084 iarg = 1 ;
00085 while( iarg < argc && argv[iarg][0] == '-' ){
00086
00087 if( strncmp(argv[iarg],"-verbose",5) == 0 ){
00088 verbose = 1 ;
00089 iarg++ ; continue ;
00090 }
00091
00092 if( strcmp(argv[iarg],"-metric") == 0 ){
00093 metric_thresh = strtod( argv[++iarg] , NULL ) ;
00094 if( metric_thresh <= 1.0 ){fprintf(stderr,"** Illegal value after -metric\n");exit(1);}
00095 iarg++ ; continue ;
00096 }
00097
00098 if( strcmp(argv[iarg],"-area") == 0 ){
00099 area_thresh = strtod( argv[++iarg] , NULL ) ;
00100 if( area_thresh <= 1.0 ){fprintf(stderr,"** Illegal value after -area\n");exit(1);}
00101 iarg++ ; continue ;
00102 }
00103
00104 if( strcmp(argv[iarg],"-base") == 0 ){
00105 base = strtol( argv[++iarg] , NULL , 10 ) ;
00106 if( base <= 0 ){fprintf(stderr,"** Illegal value after -base\n");exit(1);}
00107 iarg++ ; continue ;
00108 }
00109
00110 if( strcmp(argv[iarg],"-prefix") == 0 ){
00111 strcpy( prefix , argv[++iarg] ) ;
00112 iarg++ ; continue ;
00113 }
00114
00115 if( strcmp(argv[iarg],"-hyper") == 0 ){
00116 hyper = 1 ;
00117 iarg++ ; continue ;
00118 }
00119
00120 if( strcmp(argv[iarg],"-hyclip") == 0 ){
00121 hyper = 2 ;
00122 iarg++ ; continue ;
00123 }
00124
00125 if( strcmp(argv[iarg],"-hyval") == 0 ){
00126 hyval = strtol( argv[++iarg] , NULL , 10 ) ;
00127 if( hyval <= 0 ){fprintf(stderr,"** Illegal value after -hyval\n");exit(1);}
00128 iarg++ ; continue ;
00129 }
00130
00131 fprintf(stderr,"** Illegal option: %s\n",argv[iarg]) ; exit(1) ;
00132 }
00133
00134
00135
00136 if( iarg >= argc ){fprintf(stderr,"** No dataset name on command line?\n");exit(1);}
00137
00138 dset = THD_open_dataset( argv[iarg] ) ;
00139 if( dset == NULL ){fprintf(stderr,"** Can't open dataset %s\n",argv[iarg]);exit(1);}
00140
00141 dx = DSET_DX(dset) ; dy = DSET_DY(dset) ; dz = DSET_DZ(dset) ;
00142 nx = DSET_NX(dset) ; ny = DSET_NY(dset) ; nz = DSET_NZ(dset) ;
00143
00144 nxy = nx*ny ; nxyz = nx*ny*nz ;
00145
00146 if( dx <= 0.0 || dx != dy || dx != dz ){
00147 fprintf(stderr,"** Dataset voxel shape is illegal!\n") ; exit(1) ;
00148 }
00149
00150 if( dset->view_type == VIEW_ORIGINAL_TYPE ){
00151 fprintf(stderr,"** Dataset is in the +orig view, which is illegal!\n") ; exit(1) ;
00152 }
00153
00154 DSET_mallocize(dset) ; DSET_load(dset) ;
00155 if( !DSET_LOADED(dset) ){fprintf(stderr,"** Can't load dataset %s\n",argv[iarg]);exit(1);}
00156
00157 dar = DSET_ARRAY(dset,0) ;
00158
00159
00160
00161 if( base <= 0 || (hyval <= 0 && hyper) ){
00162 int bb , hh ;
00163 find_base_value( nxyz , dar , &bb , &hh ) ;
00164 if( base <= 0 ){
00165 base = bb ;
00166 if( verbose ) printf("-- Base value set to %d\n",base) ;
00167 }
00168 if( hyval <= 0 && hyper && hh > base+1 ){
00169 hyval = hh ;
00170 if( verbose ) printf("-- Hyperintensity value set to %d\n",hyval) ;
00171 }
00172 }
00173
00174
00175
00176 qset = EDIT_empty_copy( dset ) ;
00177 EDIT_dset_items( qset,
00178 ADN_prefix , prefix ,
00179 ADN_nvals , 1 ,
00180 ADN_ntt , 0 ,
00181 ADN_none ) ;
00182
00183
00184
00185 qar = (short *) malloc( sizeof(short) * nxyz ) ;
00186 if( qar == NULL ){fprintf(stderr,"** Can't malloc workspace!\n");exit(1);}
00187
00188 EDIT_substitute_brick( qset , 0 , MRI_short , qar ) ;
00189
00190 for( ii=pp=0 ; ii < nxyz ; ii++ ){
00191 if( dar[ii] > base ){
00192 qar[ii] = dar[ii] ;
00193 } else {
00194 qar[ii] = 0 ; pp++ ;
00195 }
00196 }
00197 if( verbose ) printf("-- Blasted %d voxels at or below base\n",pp) ;
00198
00199
00200
00201 for( jj=0 ; jj < ny ; jj++ )
00202 for( ii=0 ; ii < nx ; ii++ ){ QQ(ii,jj,0) = QQ(ii,jj,nz-1) = 0 ; }
00203
00204 for( kk=0 ; kk < nz ; kk++ )
00205 for( jj=0 ; jj < ny ; jj++ ){ QQ(0,jj,kk) = QQ(nx-1,jj,kk) = 0 ; }
00206
00207 for( kk=0 ; kk < nz ; kk++ )
00208 for( ii=0 ; ii < nx ; ii++ ){ QQ(ii,0,kk) = QQ(ii,ny-1,kk) = 0 ; }
00209
00210
00211
00212
00213 if( dset->view_type == VIEW_TALAIRACH_TYPE ){
00214 xyz_to_ijk( dset , 0.0 , 0.0 , 75.0 , NULL,NULL , &kc ) ;
00215 for( kk=kc ; kk < nz ; kk++ )
00216 for( jj=0 ; jj < ny ; jj++ )
00217 for( ii=0 ; ii < nx ; ii++ ) QQ(ii,jj,kk) = 0 ;
00218 ktop = kc-1 ;
00219 } else {
00220 ktop = nz-1 ;
00221 }
00222
00223
00224
00225 remove_isolated_stuff( qset ) ;
00226
00227
00228
00229 xyz_to_ijk( dset , 0.0,15.0,0.0 , &ic , &jc , NULL ) ;
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239 bar = (byte *) malloc( sizeof(byte) * nxy ) ;
00240 didit = 0 ;
00241
00242 for( kk=ktop ; kk >= 0 ; kk-- ){
00243
00244 memset( bar , 0 , sizeof(byte) * nxy ) ;
00245 cc = 0 ;
00246 for( jj=0 ; jj < ny ; jj++ )
00247 for( ii=0 ; ii < nx ; ii++ )
00248 if( QQ(ii,jj,kk) == 0 ){ BB(ii,jj) = 1 ; cc++ ; }
00249
00250 if( cc < 0.1*nxy ){
00251 if( verbose ) printf("-- Slice %3d: skipping because not enough zeros\n",kk) ;
00252 continue ;
00253 }
00254
00255
00256
00257 for( ii=0 ; ii < ic ; ii++ )
00258 if( BB(ii,jc) != 1 ) break ;
00259 if( ii < ic ) DRAW_2dfiller( nx,ny , ii,jc , bar ) ;
00260
00261
00262
00263 for( ii=nx-1 ; ii > ic ; ii-- )
00264 if( BB(ii,jc) != 1 ) break ;
00265 if( ii > ic && BB(ii,jc) == 0 ) DRAW_2dfiller( nx,ny , ii,jc , bar ) ;
00266
00267
00268
00269 pp = MIN(ic,jc) ;
00270 for( ii=ic-pp,jj=jc-pp ; ii < ic && jj < jc ; ii++,jj++ )
00271 if( BB(ii,jj) != 1 ) break ;
00272 if( ii < ic && jj < jc && BB(ii,jj) == 0 ) DRAW_2dfiller( nx,ny , ii,jj , bar ) ;
00273
00274 pp = MIN(ic,ny-1-jc) ;
00275 for( ii=ic-pp,jj=jc+pp ; ii < ic && jj > jc ; ii++,jj-- )
00276 if( BB(ii,jj) != 1 ) break ;
00277 if( ii < ic && jj > jc && BB(ii,jj) == 0 ) DRAW_2dfiller( nx,ny , ii,jj , bar ) ;
00278
00279 pp = MIN(nx-1-ic,jc) ;
00280 for( ii=ic+pp,jj=jc-pp ; ii > ic && jj < jc ; ii--,jj++ )
00281 if( BB(ii,jj) != 1 ) break ;
00282 if( ii > ic && jj < jc && BB(ii,jj) == 0 ) DRAW_2dfiller( nx,ny , ii,jj , bar ) ;
00283
00284 pp = MIN(nx-1-ic,ny-1-jc) ;
00285 for( ii=ic+pp,jj=jc+pp ; ii > ic && jj > jc ; ii--,jj-- )
00286 if( BB(ii,jj) != 1 ) break ;
00287 if( ii > ic && jj > jc && BB(ii,jj) == 0 ) DRAW_2dfiller( nx,ny , ii,jj , bar ) ;
00288
00289
00290
00291
00292
00293
00294
00295
00296 qsum = sum = 0.0 ; nsum = 0 ;
00297 for( jj=0 ; jj < ny ; jj++ ){
00298 for( ii=0 ; ii < nx ; ii++ ){
00299 if( BB(ii,jj) == 2 ){
00300 rr = (ii-ic)*(ii-ic) + (jj-jc)*(jj-jc) ;
00301 qsum += rr ;
00302 sum += sqrt(rr) ;
00303 nsum++ ;
00304 }
00305 }
00306 }
00307
00308 if( nsum < 10 ){
00309 if( verbose ) printf("-- Slice %3d: skipping because not enough boundary fillin\n",kk) ;
00310 continue ;
00311 }
00312
00313 qsum = qsum/nsum ;
00314 sum = (sum*sum)/(nsum*nsum) ;
00315 metric = qsum / sum ;
00316 area = nsum * dx*dx ;
00317 if( verbose ) printf("-- Slice %3d: metric =%8.4f area =%8.0f\n",kk,metric,area) ;
00318
00319 if( metric <= metric_thresh && area <= area_thresh ){
00320 for( jj=0 ; jj < ny ; jj++ ){
00321 for( ii=0 ; ii < nx ; ii++ ){
00322 if( BB(ii,jj) == 2 ) QQ(ii,jj,kk) = 0 ;
00323 }
00324 }
00325 didit++ ;
00326 } else if( didit > 0 ) {
00327 break ;
00328 }
00329 }
00330
00331 if( ! didit ){fprintf(stderr,"** No slices were trimmed -- end of run!\n");exit(1);}
00332
00333
00334
00335
00336 if( dset->view_type == VIEW_TALAIRACH_TYPE ){
00337
00338 xyz_to_ijk( dset , 0.0, -71.0, 0.0 , NULL , &jc , NULL ) ;
00339 for( kk=0 ; kk < nz ; kk++ )
00340 for( jj=0 ; jj <= jc ; jj++ )
00341 for( ii=0 ; ii < nx ; ii++ ) QQ(ii,jj,kk) = 0 ;
00342
00343 xyz_to_ijk( dset , 0.0, 103.0, 0.0 , NULL , &jc , NULL ) ;
00344 for( kk=0 ; kk < nz ; kk++ )
00345 for( jj=jc ; jj < ny ; jj++ )
00346 for( ii=0 ; ii < nx ; ii++ ) QQ(ii,jj,kk) = 0 ;
00347
00348 xyz_to_ijk( dset , -69.0, 0.0, 0.0 , &ic , NULL , NULL ) ;
00349 for( kk=0 ; kk < nz ; kk++ )
00350 for( jj=0 ; jj < ny ; jj++ )
00351 for( ii=0 ; ii <= ic ; ii++ ) QQ(ii,jj,kk) = 0 ;
00352
00353 xyz_to_ijk( dset , 69.0, 0.0, 0.0 , &ic , NULL , NULL ) ;
00354 for( kk=0 ; kk < nz ; kk++ )
00355 for( jj=0 ; jj < ny ; jj++ )
00356 for( ii=ic ; ii < nx ; ii++ ) QQ(ii,jj,kk) = 0 ;
00357
00358 xyz_to_ijk( dset , 0.0, 0.0, -43.0 , NULL , &jc , &kc ) ;
00359 for( kk=0 ; kk <= kc ; kk++ )
00360 for( jj=0 ; jj <= jc ; jj++ )
00361 for( ii=0 ; ii < nx ; ii++ ) QQ(ii,jj,kk) = 0 ;
00362 }
00363
00364
00365
00366 if( hyper && hyval > base+1 ){
00367 short val ;
00368
00369 if( hyper == 1 ) val = 0 ;
00370 else if( hyper == 2 ) val = hyval ;
00371
00372 for( ii=pp=0 ; ii < nxyz ; ii++ ){
00373 if( qar[ii] >= hyval ){ qar[ii] = val ; pp++ ; }
00374 }
00375 if( verbose ) printf("-- Processed %d voxels at or above hyval\n",pp) ;
00376 }
00377
00378
00379
00380 remove_isolated_stuff( qset ) ;
00381
00382
00383
00384 DSET_write(qset) ;
00385 exit(0) ;
00386 }
00387
00388
00389
00390 #define NPMAX 128
00391
00392 void find_base_value( int nxyz , short * sfim , int * nbase , int * hyval )
00393 {
00394 float bper = 60.0 , bmin = 1 ;
00395
00396 int ii , kk , nbin , sval , sum , nbot , a,b,c , npeak,ntop , nvox ;
00397 int * fbin ;
00398 int kmin[NPMAX] , kmax[NPMAX] ;
00399 int nmin , nmax ;
00400
00401
00402
00403 fbin = (int *) malloc( sizeof(int) * 32768 ) ;
00404 for( kk=0 ; kk < 32768 ; kk++ ) fbin[kk] = 0 ;
00405
00406 nvox = 0 ;
00407
00408 for( ii=0 ; ii < nxyz ; ii++ ){
00409 kk = sfim[ii] ; if( kk >= 0 ){ fbin[kk]++ ; nvox++ ; }
00410 }
00411
00412
00413
00414 for( kk=32767 ; kk > 0 ; kk-- ) if( fbin[kk] > 0 ) break ;
00415 if( kk == 0 ){fprintf(stderr,"** find_base_value: All voxels are zero!\n");exit(1);}
00416 nbin = kk+1 ;
00417
00418
00419
00420 sval = 0.01 * bper * nvox ;
00421 sum = 0 ;
00422 for( kk=0 ; kk < nbin ; kk++ ){
00423 sum += fbin[kk] ; if( sum >= sval ) break ;
00424 }
00425 nbot = kk ; if( nbot == 0 ) nbot = 1 ; if( bmin > nbot ) nbot = bmin ;
00426 if( nbot >= nbin-9 ){
00427 fprintf(stderr,"** find_base_value: Base point on histogram too high\n");
00428 exit(1);
00429 }
00430
00431
00432
00433 b = fbin[nbot-1] ; c = fbin[nbot] ;
00434 for( kk=nbot ; kk < nbin ; kk++ ){
00435 a = b ; b = c ; c = fbin[kk+1] ; fbin[kk] = 0.25*(a+c+2*b) ;
00436 }
00437
00438
00439
00440 nmin = nmax = 0 ;
00441 for( kk=nbot+1 ; kk < nbin ; kk++ ){
00442 if( fbin[kk] < fbin[kk-1] && fbin[kk] < fbin[kk+1] && nmin < NPMAX ){
00443 kmin[nmin++] = kk ;
00444 } else if( fbin[kk] > fbin[kk-1] && fbin[kk] > fbin[kk+1] && nmax < NPMAX ){
00445 kmax[nmax++] = kk ;
00446 }
00447 }
00448
00449
00450
00451 if( nmax == 0 ){
00452 fprintf(stderr,"** find_base_value: No histogram maxima above base point\n");
00453 exit(1);
00454 }
00455
00456 if( nmax == 1 ){
00457 npeak = kmax[0] ; ntop = 0 ;
00458 } else {
00459 int f1,f2 , k1,k2 , fk , klow,kup ;
00460
00461 k1 = 0 ; f1 = fbin[kmax[0]] ;
00462 k2 = 1 ; f2 = fbin[kmax[1]] ;
00463 if( f1 < f2 ){
00464 k1 = 1 ; f1 = fbin[kmax[1]] ;
00465 k2 = 0 ; f2 = fbin[kmax[0]] ;
00466 }
00467
00468 for( kk=2 ; kk < nmax ; kk++ ){
00469 fk = fbin[kmax[kk]] ;
00470 if( fk > f1 ){
00471 f2 = f1 ; k2 = k1 ;
00472 f1 = fk ; k1 = kk ;
00473 } else if( fk > f2 ){
00474 f2 = fk ; k2 = kk ;
00475 }
00476 }
00477 npeak = MIN( kmax[k1] , kmax[k2] ) ;
00478
00479
00480
00481 ntop = MAX( kmax[k1] , kmax[k2] ) ;
00482
00483 fk = fbin[ntop] ; klow = ntop ;
00484 for( kk=ntop-1 ; kk >= npeak ; kk-- ){
00485 if( fbin[kk] < fk ){ fk = fbin[kk] ; klow = kk ; }
00486 }
00487 fk = MAX( 0.10*fk , 0.05*fbin[ntop] ) ;
00488 kup = MIN( nbin-1 , ntop+3*(ntop-klow+2) ) ;
00489 for( kk=ntop+1 ; kk <= kup ; kk++ ) if( fbin[kk] < fk ) break ;
00490
00491 ntop = kk ;
00492 }
00493
00494 for( kk=npeak-1 ; kk > 0 ; kk-- )
00495 if( fbin[kk] < fbin[kk-1] && fbin[kk] < fbin[kk+1] ) break ;
00496
00497 if( ntop == 0 ) ntop = npeak + (npeak-kk) ;
00498
00499 *nbase = kk ;
00500 *hyval = ntop ;
00501
00502 free(fbin) ; return ;
00503 }
00504
00505
00506
00507 void remove_isolated_stuff( THD_3dim_dataset * qset )
00508 {
00509 short * qar = DSET_ARRAY(qset,0) ;
00510 short nb[27] ;
00511 int nblast , nx,ny,nz , nxy,nxyz , ii,jj,kk,ll,cc ;
00512 float dx = DSET_DX(qset) ;
00513 EDIT_options edopt ;
00514
00515 nx = DSET_NX(qset) ; ny = DSET_NY(qset) ; nz = DSET_NZ(qset) ;
00516 nxy = nx*ny ; nxyz = nx*ny*nz ;
00517
00518 do{ nblast = 0 ;
00519
00520
00521
00522 for( kk=1 ; kk < nz-1 ; kk++ ){
00523 for( jj=1 ; jj < ny-1 ; jj++ ){
00524 for( ii=1 ; ii < nx-1 ; ii++ ){
00525 if( QQ(ii,jj,kk) > 0 ){
00526 nb[ 0] = QQ(ii-1,jj-1,kk-1) ;
00527 nb[ 1] = QQ(ii ,jj-1,kk-1) ;
00528 nb[ 2] = QQ(ii+1,jj-1,kk-1) ;
00529 nb[ 3] = QQ(ii-1,jj ,kk-1) ;
00530 nb[ 4] = QQ(ii ,jj ,kk-1) ;
00531 nb[ 5] = QQ(ii+1,jj ,kk-1) ;
00532 nb[ 6] = QQ(ii-1,jj+1,kk-1) ;
00533 nb[ 7] = QQ(ii ,jj+1,kk-1) ;
00534 nb[ 8] = QQ(ii+1,jj+1,kk-1) ;
00535 nb[ 9] = QQ(ii-1,jj-1,kk ) ;
00536 nb[10] = QQ(ii ,jj-1,kk ) ;
00537 nb[11] = QQ(ii+1,jj-1,kk ) ;
00538 nb[12] = QQ(ii-1,jj ,kk ) ;
00539 nb[13] = QQ(ii ,jj ,kk ) ;
00540 nb[14] = QQ(ii+1,jj ,kk ) ;
00541 nb[15] = QQ(ii-1,jj-1,kk ) ;
00542 nb[16] = QQ(ii ,jj-1,kk ) ;
00543 nb[17] = QQ(ii+1,jj-1,kk ) ;
00544 nb[18] = QQ(ii-1,jj+1,kk+1) ;
00545 nb[19] = QQ(ii ,jj+1,kk+1) ;
00546 nb[20] = QQ(ii+1,jj+1,kk+1) ;
00547 nb[21] = QQ(ii-1,jj ,kk+1) ;
00548 nb[22] = QQ(ii ,jj ,kk+1) ;
00549 nb[23] = QQ(ii+1,jj ,kk+1) ;
00550 nb[24] = QQ(ii-1,jj+1,kk+1) ;
00551 nb[25] = QQ(ii ,jj+1,kk+1) ;
00552 nb[26] = QQ(ii+1,jj+1,kk+1) ;
00553
00554 for( ll=cc=0 ; ll < 27 ; ll++ ) if( nb[ll] > 0 ) cc++ ;
00555 if( cc < 4 ){ QQ(ii,jj,kk) = 0 ; nblast++ ; }
00556 }
00557 }
00558 }
00559 }
00560
00561
00562
00563 for( kk=1 ; kk < nz-1 ; kk++ ){
00564 for( jj=1 ; jj < ny-1 ; jj++ ){
00565 for( ii=1 ; ii < nx-1 ; ii++ ){
00566 if( QQ(ii,jj,kk) > 0 ){
00567 nb[ 0] = QQ(ii-1,jj-1,kk) ;
00568 nb[ 1] = QQ(ii ,jj-1,kk) ;
00569 nb[ 2] = QQ(ii+1,jj-1,kk) ;
00570 nb[ 3] = QQ(ii-1,jj ,kk) ;
00571 nb[ 4] = QQ(ii ,jj ,kk) ;
00572 nb[ 5] = QQ(ii+1,jj ,kk) ;
00573 nb[ 6] = QQ(ii-1,jj+1,kk) ;
00574 nb[ 7] = QQ(ii ,jj+1,kk) ;
00575 nb[ 8] = QQ(ii+1,jj+1,kk) ;
00576 for( ll=cc=0 ; ll < 9 ; ll++ ) if( nb[ll] > 0 ) cc++ ;
00577 if( cc < 2 ){ QQ(ii,jj,kk) = 0 ; nblast++ ; }
00578 }
00579 }
00580 }
00581 }
00582
00583 for( jj=1 ; jj < ny-1 ; jj++ ){
00584 for( kk=1 ; kk < nz-1 ; kk++ ){
00585 for( ii=1 ; ii < nx-1 ; ii++ ){
00586 if( QQ(ii,jj,kk) > 0 ){
00587 nb[ 0] = QQ(ii-1,jj,kk-1) ;
00588 nb[ 1] = QQ(ii ,jj,kk-1) ;
00589 nb[ 2] = QQ(ii+1,jj,kk-1) ;
00590 nb[ 3] = QQ(ii-1,jj,kk ) ;
00591 nb[ 4] = QQ(ii ,jj,kk ) ;
00592 nb[ 5] = QQ(ii+1,jj,kk ) ;
00593 nb[ 6] = QQ(ii-1,jj,kk+1) ;
00594 nb[ 7] = QQ(ii ,jj,kk+1) ;
00595 nb[ 8] = QQ(ii+1,jj,kk+1) ;
00596 for( ll=cc=0 ; ll < 9 ; ll++ ) if( nb[ll] > 0 ) cc++ ;
00597 if( cc < 2 ){ QQ(ii,jj,kk) = 0 ; nblast++ ; }
00598 }
00599 }
00600 }
00601 }
00602
00603 for( ii=1 ; ii < nx-1 ; ii++ ){
00604 for( jj=1 ; jj < ny-1 ; jj++ ){
00605 for( kk=1 ; kk < nz-1 ; kk++ ){
00606 if( QQ(ii,jj,kk) > 0 ){
00607 nb[ 0] = QQ(ii,jj-1,kk-1) ;
00608 nb[ 1] = QQ(ii,jj ,kk-1) ;
00609 nb[ 2] = QQ(ii,jj+1,kk-1) ;
00610 nb[ 3] = QQ(ii,jj-1,kk ) ;
00611 nb[ 4] = QQ(ii,jj ,kk ) ;
00612 nb[ 5] = QQ(ii,jj+1,kk ) ;
00613 nb[ 6] = QQ(ii,jj-1,kk+1) ;
00614 nb[ 7] = QQ(ii,jj ,kk+1) ;
00615 nb[ 8] = QQ(ii,jj+1,kk+1) ;
00616 for( ll=cc=0 ; ll < 9 ; ll++ ) if( nb[ll] > 0 ) cc++ ;
00617 if( cc < 2 ){ QQ(ii,jj,kk) = 0 ; nblast++ ; }
00618 }
00619 }
00620 }
00621 }
00622
00623 if( verbose ) printf("-- Blasted %d voxels for being too isolated\n",nblast) ;
00624 } while( nblast > 0 ) ;
00625
00626
00627
00628 INIT_EDOPT( &edopt ) ;
00629 edopt.edit_clust = ECFLAG_SAME;
00630 edopt.clust_rmm = 1.01*dx ;
00631 edopt.clust_vmul = 25000.0 ;
00632 EDIT_one_dataset( qset , &edopt ) ;
00633
00634 return ;
00635 }
00636
00637
00638
00639 void xyz_to_ijk( THD_3dim_dataset * ds , float x , float y , float z ,
00640 int * i , int * j , int * k )
00641 {
00642 THD_fvec3 fv ;
00643 THD_ivec3 iv ;
00644
00645 LOAD_FVEC3( fv , x,y,z ) ;
00646 fv = THD_dicomm_to_3dmm( ds , fv ) ;
00647 iv = THD_3dmm_to_3dind ( ds , fv ) ;
00648 if( i != NULL ) *i = iv.ijk[0] ;
00649 if( j != NULL ) *j = iv.ijk[1] ;
00650 if( k != NULL ) *k = iv.ijk[2] ;
00651 return ;
00652 }
00653
00654
00655
00656
00657
00658
00659
00660
00661
00662
00663
00664
00665 void DRAW_2dfiller( int nx , int ny , int ix , int jy , byte * ar )
00666 {
00667 int ii,jj , ip,jp , num ;
00668
00669 #define AR(i,j) ar[(i)+(j)*nx]
00670
00671
00672
00673 ip = ix ; jp = jy ; AR(ip,jp) = 2 ;
00674
00675 for( ii=ip+1; ii < nx && AR(ii,jp) == 0; ii++ ) AR(ii,jp) = 2;
00676 for( ii=ip-1; ii >= 0 && AR(ii,jp) == 0; ii-- ) AR(ii,jp) = 2;
00677 for( jj=jp+1; jj < ny && AR(ip,jj) == 0; jj++ ) AR(ip,jj) = 2;
00678 for( jj=jp-1; jj >= 0 && AR(ip,jj) == 0; jj-- ) AR(ip,jj) = 2;
00679
00680
00681
00682 do {
00683 num = 0 ;
00684 for( jp=0 ; jp < ny ; jp++ ){
00685 for( ip=0 ; ip < nx ; ip++ ){
00686 if( AR(ip,jp) == 2 ){
00687 for( ii=ip+1; ii < nx && AR(ii,jp) == 0; ii++ ){ AR(ii,jp) = 2; num++; }
00688 for( ii=ip-1; ii >= 0 && AR(ii,jp) == 0; ii-- ){ AR(ii,jp) = 2; num++; }
00689 for( jj=jp+1; jj < ny && AR(ip,jj) == 0; jj++ ){ AR(ip,jj) = 2; num++; }
00690 for( jj=jp-1; jj >= 0 && AR(ip,jj) == 0; jj-- ){ AR(ip,jj) = 2; num++; }
00691 }
00692 }
00693 }
00694 } while( num > 0 ) ;
00695
00696 return ;
00697 }