Doxygen Source Code Documentation
3dstrip.c File Reference
#include "mrilib.h"Go to the source code of this file.
Defines | |
| #define | DD(i, j, k) dar[(i)+(j)*nx+(k)*nxy] |
| #define | QQ(i, j, k) qar[(i)+(j)*nx+(k)*nxy] |
| #define | BB(i, j) bar[(i)+(j)*nx] |
| #define | NPMAX 128 |
| #define | AR(i, j) ar[(i)+(j)*nx] |
Functions | |
| void | find_base_value (int nxyz, short *sfim, int *base, int *peak) |
| void | remove_isolated_stuff (THD_3dim_dataset *qset) |
| void | xyz_to_ijk (THD_3dim_dataset *ds, float x, float y, float z, int *i, int *j, int *k) |
| void | DRAW_2dfiller (int nx, int ny, int ix, int jy, byte *ar) |
| int | main (int argc, char *argv[]) |
Variables | |
| int | verbose = 0 |
Define Documentation
|
|
|
|
|
Definition at line 19 of file 3dstrip.c. Referenced by main(). |
|
|
|
|
|
Definition at line 390 of file 3dstrip.c. Referenced by find_base_value(). |
|
|
Definition at line 18 of file 3dstrip.c. Referenced by main(), and remove_isolated_stuff(). |
Function Documentation
|
||||||||||||||||||||||||
|
Definition at line 665 of file 3dstrip.c.
00666 {
00667 int ii,jj , ip,jp , num ;
00668
00669 #define AR(i,j) ar[(i)+(j)*nx]
00670
00671 /* fill out in cross from 1st point */
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 /* brute force repetition of the cross technique */
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 }
|
|
||||||||||||||||||||
|
Definition at line 392 of file 3dstrip.c. References a, c, free, malloc, MAX, MIN, and NPMAX.
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 /*-- make histogram of shorts --*/
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 /*-- find largest value --*/
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 /*-- find bper point in cumulative distribution --*/
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 /*-- smooth histogram --*/
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 /*-- find minima and maxima above bper point --*/
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 /*-- find the largest two maxima --*/
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] ) ; /* smaller bin of the 2 top peaks */
00478
00479 /* find valley between 2 peaks */
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 }
|
|
||||||||||||
|
---------- Adapted from 3dZeropad.c by RWCox - 08 Aug 2001 ----------* Definition at line 23 of file 3dstrip.c. References ADN_none, ADN_ntt, ADN_nvals, ADN_prefix, argc, base, BB, DRAW_2dfiller(), DSET_ARRAY, DSET_DX, DSET_DY, DSET_DZ, DSET_load, DSET_LOADED, DSET_mallocize, DSET_NX, DSET_NY, DSET_NZ, DSET_write, EDIT_dset_items(), EDIT_empty_copy(), EDIT_substitute_brick(), find_base_value(), malloc, MIN, nz, QQ, remove_isolated_stuff(), strtod(), THD_open_dataset(), verbose, VIEW_ORIGINAL_TYPE, VIEW_TALAIRACH_TYPE, THD_3dim_dataset::view_type, and xyz_to_ijk().
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 /* help? */
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 /* arguments */
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 /* open dataset */
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 /* find base value, if needed */
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 /* make new dataset */
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 /* copy dar to qar, omiting values <= base */
00184
00185 qar = (short *) malloc( sizeof(short) * nxyz ) ; /* another array */
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 /* set edges to zero */
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 /* if in Talairach view, remove stuff above the brain */
00211 /* [coordinate order is x=R-L, y=A-P, z=I-S] */
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 /* now prune the dataset of the little isolated junk */
00224
00225 remove_isolated_stuff( qset ) ;
00226
00227 /* find center of brain */
00228
00229 xyz_to_ijk( dset , 0.0,15.0,0.0 , &ic , &jc , NULL ) ; /* a decent place */
00230
00231 /* for each slice:
00232 build a mask of blocking points (zero voxels in the slice)
00233 find some outermost points
00234 flood fill from the outermost points
00235 compute a metric of "annularity" of the filled points
00236 if the slice passes the metric, blast out the filled points
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 ) ; /* build the mask */
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 /* scan in from right, fill from the first nonmasked point we hit */
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 /* scan in from left */
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 /* scan in from diagonals */
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 /* we don't scan in from top/bottom since those might be
00290 clipped off in Talairach coordinates, which would be disastrous! */
00291
00292 /* compute metric for marked points:
00293 metric = average( r**2 ) / average(r)**2
00294 this should be nearly 1 for an annular region */
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 ; /* don't go below first slice that fails, after some slice worked */
00328 }
00329 }
00330
00331 if( ! didit ){fprintf(stderr,"** No slices were trimmed -- end of run!\n");exit(1);}
00332
00333 /* if in Talairach view, remove stuff to the sides of the brain */
00334 /* (done after the slice stuff so the scalp annulus isn't broken) */
00335
00336 if( dset->view_type == VIEW_TALAIRACH_TYPE ){
00337
00338 xyz_to_ijk( dset , 0.0, -71.0, 0.0 , NULL , &jc , NULL ) ; /* anterior */
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 ) ; /* posterior */
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 ) ; /* right */
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 ) ; /* left */
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 ) ; /* below */
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 /* get rid of the hyperintensity stuff */
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 /* re-prune the dataset of the little isolated junk */
00379
00380 remove_isolated_stuff( qset ) ;
00381
00382 /* write output */
00383
00384 DSET_write(qset) ;
00385 exit(0) ;
00386 }
|
|
|
Definition at line 507 of file 3dstrip.c. References EDIT_options::clust_rmm, EDIT_options::clust_vmul, DSET_ARRAY, DSET_DX, DSET_NX, DSET_NY, DSET_NZ, ECFLAG_SAME, EDIT_options::edit_clust, EDIT_one_dataset(), INIT_EDOPT, nz, and QQ. Referenced by main().
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 /* edit 3x3x3 volumes */
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 ){ /* must have at least 3 nonzero */
00526 nb[ 0] = QQ(ii-1,jj-1,kk-1) ; /* voxels in a 3x3x3 neighborhood */
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 /* edit 3x3 areas in each slice orientation */
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 ){ /* must have at least 2 nonzero */
00567 nb[ 0] = QQ(ii-1,jj-1,kk) ; /* voxels in a 3x3 neighborhood */
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 ){ /* must have at least 2 nonzero */
00587 nb[ 0] = QQ(ii-1,jj,kk-1) ; /* voxels in a 3x3 neighborhood */
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 ){ /* must have at least 2 nonzero */
00607 nb[ 0] = QQ(ii,jj-1,kk-1) ; /* voxels in a 3x3 neighborhood */
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 /* cluster data, remove small clusters */
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 }
|
|
||||||||||||||||||||||||||||||||
|
Definition at line 639 of file 3dstrip.c. References i, THD_ivec3::ijk, LOAD_FVEC3, THD_3dmm_to_3dind(), and THD_dicomm_to_3dmm(). Referenced by main().
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 }
|
Variable Documentation
|
|
Definition at line 21 of file 3dstrip.c. Referenced by main(). |
