Doxygen Source Code Documentation
edt_onedset.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Defines | |
#define | THADD |
Functions | |
void | EDIT_one_dataset (THD_3dim_dataset *dset, EDIT_options *edopt) |
Define Documentation
|
|
Function Documentation
|
Definition at line 22 of file edt_onedset.c. References EDIT_options::abss, ADDTO_CLARR, ADN_brick_stataux_one, ADN_none, AFNI_GOOD_DTYPE, AFNI_GOOD_FUNC_DTYPE, EDIT_options::blur, THD_diskptr::brick_name, CABS, MCW_cluster_array::clar, EDIT_options::clip_bot, EDIT_options::clip_top, EDIT_options::clip_unscaled, EDIT_options::clust_rmm, EDIT_options::clust_vmul, DATABLOCK_MEM_MALLOC, THD_3dim_dataset::daxes, THD_3dim_dataset::dblk, DESTROY_CLARR, EDIT_options::dilate, THD_datablock::diskptr, EDIT_options::do_zvol, DSET_ARRAY, DSET_BRICK_FACTOR, DSET_BRICK_STATAUX, DSET_BRICK_STATCODE, DSET_BRICK_TYPE, DSET_BRIKNAME, DSET_LOADED, DSET_NVALS, ECFLAG_ORDER, ECFLAG_SAME, ECFLAG_SIZE, EDIT_blur_volume(), EDIT_options::edit_clust, EDIT_cluster_array(), EDIT_coerce_scale_type(), EDIT_dset_items(), EDIT_filter_volume(), EDIT_zscore_vol(), ENTRY, EDIT_options::erode_pv, EDIT_options::fake_dxyz, FCFLAG_NONE, EDIT_options::fexpr, EDIT_options::filter_opt, EDIT_options::filter_rmm, EDIT_options::fmask, FUNC_IS_STAT, THD_3dim_dataset::func_type, FUNC_ZT_SCALE_SHORT, complex::i, THD_ivec3::ijk, INIT_CLARR, ISANAT, ISBUCKET, ISFUNC, EDIT_options::iv_fim, EDIT_options::iv_thr, MAX, MCW_cluster_to_vol(), MCW_erode_clusters(), MCW_find_clusters(), MCW_scale_to_max(), MRI_type_name, EDIT_options::mult, EDIT_options::noneg, MCW_cluster_array::num_clu, MCW_cluster::num_pt, THD_dataxes::nxx, THD_dataxes::nyy, nz, THD_dataxes::nzz, complex::r, EDIT_options::scale, SORT_CLARR, THD_3dim_dataset::stat_aux, STATUS, TEMP_FVEC3, THD_3dmm_to_3dind(), THD_force_malloc_type(), THD_load_datablock(), EDIT_options::thrblur, EDIT_options::thresh, EDIT_options::thrfilter_opt, EDIT_options::thrfilter_rmm, EDIT_options::thtoin, top, EDIT_options::verbose, THD_dataxes::xxdel, THD_dataxes::yydel, EDIT_options::zscore, EDIT_options::zv_x1, EDIT_options::zv_x2, EDIT_options::zv_y1, EDIT_options::zv_y2, EDIT_options::zv_z1, EDIT_options::zv_z2, and THD_dataxes::zzdel. Referenced by CLUST_main(), EDIT_main(), main(), and remove_isolated_stuff().
00023 { 00024 int edit_thtoin = edopt->thtoin ; /* copy into local variables */ 00025 int edit_noneg = edopt->noneg ; /* for historical reasons */ 00026 int edit_abs = edopt->abss ; 00027 float edit_clip_bot = edopt->clip_bot ; /* Nov 1995: changed to floats */ 00028 float edit_clip_top = edopt->clip_top ; 00029 float edit_thresh = edopt->thresh ; 00030 int edit_clust = edopt->edit_clust ; /* 10 Sept 1996 */ 00031 float clust_rmm = edopt->clust_rmm ; 00032 float clust_vmul = edopt->clust_vmul ; 00033 float erode_pv = edopt->erode_pv; /* 17 June 1998 */ 00034 int dilate = edopt->dilate; /* 17 June 1998 */ 00035 int filter_opt = edopt->filter_opt; /* 11 Sept 1996 */ 00036 float filter_rmm = edopt->filter_rmm; /* 11 Sept 1996 */ 00037 int thrfilter_opt = edopt->thrfilter_opt; /* 1 Oct 1996 */ 00038 float thrfilter_rmm = edopt->thrfilter_rmm; /* 1 Oct 1996 */ 00039 float edit_blur = edopt->blur ; 00040 float edit_thrblur = edopt->thrblur; /* 4 Oct 1996 */ 00041 int edit_scale = edopt->scale ; 00042 float edit_mult = edopt->mult ; 00043 int edit_zvol = edopt->do_zvol ; 00044 int edit_ivfim = edopt->iv_fim ; /* 30 Nov 1997 */ 00045 int edit_ivthr = edopt->iv_thr ; /* 30 Nov 1997 */ 00046 int verbose = edopt->verbose ; /* 01 Nov 1999 */ 00047 int fake_dxyz = edopt->fake_dxyz ; /* 11 Sep 2000 */ 00048 00049 int edit_clip_unscaled = edopt->clip_unscaled ; /* 09 Aug 1996 */ 00050 00051 THD_dataxes * daxes ; 00052 short * sfim = NULL , * sthr = NULL ; 00053 float * ffim = NULL , * fthr = NULL ; 00054 complex * cfim = NULL ; 00055 byte * bfim = NULL , * bthr = NULL ; 00056 void * vfim = NULL , * vthr = NULL ; 00057 int nx,ny,nz,nxy,nxyz , jj,kk , ptmin , iclu,nclu , fim_max ; 00058 int iv_fim , iv_thr , fim_type , thr_type ; 00059 register int ii ; 00060 float dx,dy,dz , dxyz , rmm,vmul , val , vvv ; 00061 MCW_cluster_array * clar ; 00062 MCW_cluster * blur=NULL ; 00063 int fimtype , thrtype ; 00064 float fimfac , thrfac ; 00065 00066 /** get the data from this dataset **/ 00067 00068 ENTRY("EDIT_one_dataset") ; 00069 00070 THD_force_malloc_type( dset->dblk , DATABLOCK_MEM_MALLOC ) ; 00071 THD_load_datablock( dset->dblk ) ; 00072 00073 if( !DSET_LOADED(dset) ){ 00074 fprintf(stderr, 00075 "\n*** Cannot read data brick for dataset %s\a\n", 00076 DSET_BRIKNAME(dset) ) ; 00077 EXRETURN ; 00078 } 00079 00080 /** load the data sub-brick indexes (iv_*) and check types for legality **/ 00081 00082 if( ISANAT(dset) ){ 00083 if( edit_ivfim >= 0 && edit_ivfim < DSET_NVALS(dset) ) /* 30 Nov 1997 */ 00084 iv_fim = edit_ivfim ; 00085 else 00086 iv_fim = ANAT_ival_zero[dset->func_type] ; 00087 00088 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: fim index = %d\n",iv_fim) ; 00089 00090 fim_type = DSET_BRICK_TYPE(dset,iv_fim) ; 00091 fimfac = DSET_BRICK_FACTOR(dset,iv_fim) ; 00092 iv_thr = -1 ; 00093 thr_type = ILLEGAL_TYPE ; 00094 00095 if( !AFNI_GOOD_DTYPE(fim_type) || fim_type == MRI_rgb ){ 00096 fprintf(stderr,"\n*** Illegal anatomy data type %s in dataset %s\a\n" , 00097 MRI_type_name[fim_type] , 00098 dset->dblk->diskptr->brick_name ) ; 00099 EXRETURN ; 00100 } 00101 00102 #ifdef AFNI_DEBUG 00103 { char str[256] ; 00104 sprintf(str,"Anat dset: iv=%d type=%s fac=%g",iv_fim,MRI_TYPE_name[fim_type],fimfac) ; 00105 STATUS(str) ; } 00106 #endif 00107 00108 } 00109 00110 if( ISFUNC(dset) ){ 00111 if( edit_ivfim >= 0 && edit_ivfim < DSET_NVALS(dset) ) /* 30 Nov 1997 */ 00112 iv_fim = edit_ivfim ; 00113 else 00114 iv_fim = FUNC_ival_fim[dset->func_type] ; 00115 00116 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: fim index = %d\n",iv_fim) ; 00117 00118 fim_type = DSET_BRICK_TYPE(dset,iv_fim) ; 00119 fimfac = DSET_BRICK_FACTOR(dset,iv_fim) ; 00120 00121 if( edit_ivthr >= 0 && edit_ivthr < DSET_NVALS(dset) ) /* 30 Nov 1997 */ 00122 iv_thr = edit_ivthr ; 00123 else 00124 iv_thr = FUNC_ival_thr[dset->func_type] ; 00125 00126 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: thr index = %d\n",iv_thr) ; 00127 00128 if( iv_thr < 0 ){ 00129 thr_type = ILLEGAL_TYPE ; 00130 thrfac = 0.0 ; 00131 } else { 00132 thr_type = DSET_BRICK_TYPE(dset,iv_thr) ; 00133 thrfac = DSET_BRICK_FACTOR(dset,iv_thr) ; 00134 if( thrfac == 0.0 ){ 00135 switch( thr_type ){ 00136 case MRI_short: thrfac = 1.0/FUNC_scale_short[dset->func_type]; break; 00137 case MRI_byte : thrfac = 1.0/FUNC_scale_byte [dset->func_type]; break; 00138 } 00139 } 00140 } 00141 00142 if( !AFNI_GOOD_FUNC_DTYPE(fim_type) || fim_type == MRI_rgb ){ 00143 fprintf(stderr,"\n*** Illegal functional data type %s in dataset %s\a\n" , 00144 MRI_type_name[fim_type], dset->dblk->diskptr->brick_name ) ; 00145 EXRETURN ; 00146 } 00147 00148 if( thr_type >= 0 && (!AFNI_GOOD_FUNC_DTYPE(thr_type) || fim_type == MRI_rgb) ){ 00149 fprintf(stderr,"\n*** Illegal threshold data type %s in dataset %s\a\n" , 00150 MRI_type_name[fim_type] , dset->dblk->diskptr->brick_name ) ; 00151 EXRETURN ; 00152 } 00153 00154 #ifdef AFNI_DEBUG 00155 { char str[256] ; 00156 sprintf(str,"Func dset: iv_fim=%d type=%s fac=%g",iv_fim,MRI_TYPE_name[fim_type],fimfac) ; 00157 STATUS(str) ; 00158 if( iv_thr >= 0 ){ 00159 sprintf(str,"Func dset: iv_thr=%d type=%s fac=%g",iv_thr,MRI_TYPE_name[thr_type],thrfac) ; 00160 STATUS(str) ; } } 00161 #endif 00162 00163 } 00164 00165 /** load the pointers to the sub-bricks **/ 00166 00167 vfim = DSET_ARRAY(dset,iv_fim) ; 00168 switch( fim_type ){ 00169 default: 00170 fprintf(stderr,"\n*** Illegal data type in dataset %s\a\n", 00171 dset->dblk->diskptr->brick_name ) ; 00172 EXRETURN ; 00173 00174 case MRI_short: sfim = (short *) vfim ; break ; 00175 case MRI_float: ffim = (float *) vfim ; break ; 00176 case MRI_byte: bfim = (byte *) vfim ; break ; 00177 case MRI_complex: cfim = (complex *) vfim ; break ; 00178 } 00179 00180 if( iv_thr >= 0 ){ 00181 vthr = DSET_ARRAY(dset,iv_thr) ; 00182 switch( thr_type ){ 00183 default: 00184 fprintf(stderr,"\n*** Illegal thresh data type in dataset %s\a\n", 00185 dset->dblk->diskptr->brick_name ) ; 00186 EXRETURN ; 00187 00188 case MRI_short: sthr = (short *) vthr ; break ; 00189 case MRI_float: fthr = (float *) vthr ; break ; 00190 case MRI_byte: bthr = (byte *) vthr ; break ; 00191 } 00192 } 00193 00194 /** load the grid parameters **/ 00195 00196 daxes = dset->daxes ; 00197 nx = daxes->nxx ; dx = fabs(daxes->xxdel) ; 00198 ny = daxes->nyy ; dy = fabs(daxes->yydel) ; 00199 nz = daxes->nzz ; dz = fabs(daxes->zzdel) ; 00200 00201 if( fake_dxyz ) dx = dy = dz = 1.0 ; /* 11 Sep 2000 */ 00202 00203 nxy = nx * ny ; nxyz = nxy * nz ; dxyz = dx*dy*dz ; 00204 00205 /*----- copy threshold over intensity? -----*/ 00206 00207 STATUS("dataset loaded") ; 00208 00209 if( edit_thtoin && iv_thr >= 0 ){ 00210 float new_fimfac , scaling ; 00211 00212 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: copy thr over fim\n") ; 00213 00214 /**** 00215 Find scaling factors for various conversions (0 --> no scaling) 00216 scaling = factor to actually scale data by when copying to new brick 00217 new_fimfac = factor to later scale data by when converting to floats 00218 ****/ 00219 00220 if( edit_thtoin == 2 ){ 00221 new_fimfac = scaling = 0.0 ; /** -2thtoin --> no scaling **/ 00222 } else { 00223 switch( thr_type ){ 00224 00225 /** threshold datum is shorts **/ 00226 00227 case MRI_short:{ 00228 switch( fim_type ){ 00229 case MRI_short: /* fim datum is shorts --> no new scaling needed */ 00230 new_fimfac = thrfac ; 00231 scaling = 0.0 ; 00232 break ; 00233 00234 case MRI_float: /* fim datum is floats --> will be scaled properly */ 00235 new_fimfac = 0.0 ; 00236 scaling = thrfac ; 00237 break ; 00238 00239 case MRI_byte: /* fim datum is bytes */ 00240 new_fimfac = 1.0 / FUNC_scale_byte[dset->func_type] ; 00241 scaling = thrfac * FUNC_scale_byte[dset->func_type] ; 00242 break ; 00243 } 00244 } 00245 break ; 00246 00247 /** threshold datum is bytes **/ 00248 00249 case MRI_byte:{ 00250 switch( fim_type ){ 00251 case MRI_short: /* fim datum is shorts */ 00252 new_fimfac = 1.0 / FUNC_scale_short[dset->func_type] ; 00253 scaling = thrfac * FUNC_scale_short[dset->func_type] ; 00254 break ; 00255 00256 case MRI_float: /* fim datum is floats */ 00257 new_fimfac = 0.0 ; 00258 scaling = thrfac ; 00259 break ; 00260 00261 case MRI_byte: /* fim datum is bytes */ 00262 new_fimfac = thrfac ; 00263 scaling = 0.0 ; 00264 break ; 00265 } 00266 } 00267 break ; 00268 00269 /** threshold datum is floats **/ 00270 00271 case MRI_float:{ 00272 switch( fim_type ){ 00273 case MRI_short: /* fim datum is shorts */ 00274 new_fimfac = 1.0 / FUNC_scale_short[dset->func_type] ; 00275 scaling = FUNC_scale_short[dset->func_type] ; 00276 break ; 00277 00278 case MRI_float: /* fim datum is floats --> no scaling needed */ 00279 new_fimfac = 0.0 ; 00280 scaling = 0.0 ; 00281 break ; 00282 00283 case MRI_byte: /* fim datum is bytes */ 00284 new_fimfac = 1.0 / FUNC_scale_byte[dset->func_type] ; 00285 scaling = FUNC_scale_byte[dset->func_type] ; 00286 break ; 00287 } 00288 } 00289 break ; 00290 } 00291 } 00292 00293 #ifdef AFNI_DEBUG 00294 { char str[256] ; 00295 sprintf(str,"thtoin: scaling=%f new_fimfac=%f input=%s output=%s", 00296 scaling,new_fimfac,MRI_TYPE_name[thr_type],MRI_TYPE_name[fim_type]) ; 00297 STATUS(str) ; } 00298 #endif 00299 00300 /** have scaling factors, so use them **/ 00301 00302 EDIT_coerce_scale_type( nxyz , scaling , 00303 thr_type , vthr , fim_type , vfim ) ; 00304 00305 DSET_BRICK_FACTOR(dset,iv_fim) = fimfac = new_fimfac ; 00306 } /* end -1thtoin */ 00307 00308 /*----- non-negative? -----*/ 00309 00310 if( edit_noneg ){ /* meaningless for byte and complex */ 00311 STATUS("noneg") ; 00312 00313 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: remove negative values\n") ; 00314 00315 switch( fim_type ){ 00316 case MRI_short: 00317 for( ii=0 ; ii < nxyz ; ii++ ) if( sfim[ii] < 0 ) sfim[ii] = 0 ; 00318 break ; 00319 00320 case MRI_float: 00321 for( ii=0 ; ii < nxyz ; ii++ ) if( ffim[ii] < 0 ) ffim[ii] = 0 ; 00322 break ; 00323 00324 default: 00325 STATUS("noneg applied to meaningless type: will be ignored") ; 00326 } 00327 } 00328 00329 /*----- absolute? -----*/ 00330 00331 if( edit_abs ){ /* meaningless for byte */ 00332 STATUS("abs") ; 00333 00334 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: take absolute value\n") ; 00335 00336 switch( fim_type ){ 00337 case MRI_short: 00338 for( ii=0 ; ii < nxyz ; ii++ ) if( sfim[ii] < 0 ) sfim[ii] = -sfim[ii] ; 00339 break ; 00340 00341 case MRI_float: 00342 for( ii=0 ; ii < nxyz ; ii++ ) if( ffim[ii] < 0 ) ffim[ii] = -ffim[ii] ; 00343 break ; 00344 00345 case MRI_complex: 00346 for( ii=0 ; ii < nxyz ; ii++ ){ 00347 cfim[ii].r = CABS(cfim[ii]) ; cfim[ii].i = 0.0 ; 00348 } 00349 break ; 00350 00351 default: 00352 STATUS("abs applied to meaningless type: will be ignored") ; 00353 } 00354 } 00355 00356 /*----- clip? -----*/ 00357 00358 if( edit_clip_bot < edit_clip_top ){ 00359 00360 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: clip fim values\n") ; 00361 00362 switch( fim_type ){ 00363 case MRI_short:{ 00364 int top , bot ; 00365 float ftop,fbot ; 00366 if( fimfac > 0.0 && ! edit_clip_unscaled ){ 00367 ftop = edit_clip_top / fimfac ; 00368 fbot = edit_clip_bot / fimfac ; 00369 } else { 00370 ftop = edit_clip_top ; 00371 fbot = edit_clip_bot ; 00372 } 00373 00374 top = rint(ftop) ; /* this code was modifed 28 Sep 1998 */ 00375 if( top >= MRI_maxshort ) top = MRI_maxshort + 1 ; 00376 if( top <= -MRI_maxshort ) top = -(MRI_maxshort + 1) ; 00377 00378 bot = rint(fbot) ; 00379 if( bot >= MRI_maxshort ) bot = MRI_maxshort + 1 ; 00380 if( bot <= -MRI_maxshort ) bot = -(MRI_maxshort + 1) ; 00381 00382 #ifdef AFNI_DEBUG 00383 { char str[256] ; 00384 sprintf(str,"clipping short from %d to %d",bot,top) ; 00385 STATUS(str) ; } 00386 #endif 00387 for( ii=0 ; ii < nxyz ; ii++ ) 00388 if( sfim[ii] > bot && sfim[ii] < top ) sfim[ii] = 0 ; 00389 } 00390 break ; 00391 00392 case MRI_byte:{ 00393 int top , bot ; 00394 float ftop,fbot ; 00395 if( fimfac > 0.0 && ! edit_clip_unscaled ){ 00396 ftop = edit_clip_top / fimfac ; 00397 fbot = edit_clip_bot / fimfac ; 00398 } else { 00399 ftop = edit_clip_top ; 00400 fbot = edit_clip_bot ; 00401 } 00402 00403 top = rint(ftop) ; 00404 if( top >= MRI_maxbyte ) top = MRI_maxbyte + 1 ; 00405 if( top <= -MRI_maxbyte ) top = -(MRI_maxbyte + 1) ; 00406 00407 bot = rint(fbot) ; 00408 if( bot >= MRI_maxbyte ) bot = MRI_maxbyte + 1 ; 00409 if( bot <= -MRI_maxbyte ) bot = -(MRI_maxbyte + 1) ; 00410 00411 if( bot < 0 ) bot = 0 ; 00412 if( top < bot ) top = bot ; 00413 #ifdef AFNI_DEBUG 00414 { char str[256] ; 00415 sprintf(str,"clipping byte from %d to %d",bot,top) ; 00416 STATUS(str) ; } 00417 #endif 00418 for( ii=0 ; ii < nxyz ; ii++ ) 00419 if( bfim[ii] > bot && bfim[ii] < top ) bfim[ii] = 0 ; 00420 } 00421 break ; 00422 00423 case MRI_float:{ 00424 float top , bot ; 00425 if( fimfac > 0.0 && ! edit_clip_unscaled ){ 00426 top = edit_clip_top / fimfac ; 00427 bot = edit_clip_bot / fimfac ; 00428 } else { 00429 top = edit_clip_top ; 00430 bot = edit_clip_bot ; 00431 } 00432 #ifdef AFNI_DEBUG 00433 { char str[256] ; 00434 sprintf(str,"clipping float from %g to %g",bot,top) ; 00435 STATUS(str) ; } 00436 #endif 00437 for( ii=0 ; ii < nxyz ; ii++ ) 00438 if( ffim[ii] > bot && ffim[ii] < top ) ffim[ii] = 0.0 ; 00439 } 00440 break ; 00441 00442 case MRI_complex:{ 00443 float val ; 00444 float top , bot ; 00445 if( fimfac > 0.0 && ! edit_clip_unscaled ){ 00446 top = edit_clip_top / fimfac ; 00447 bot = edit_clip_bot / fimfac ; 00448 } else { 00449 top = edit_clip_top ; 00450 bot = edit_clip_bot ; 00451 } 00452 #ifdef AFNI_DEBUG 00453 { char str[256] ; 00454 sprintf(str,"clipping complex from %g to %g",bot,top) ; 00455 STATUS(str) ; } 00456 #endif 00457 for( ii=0 ; ii < nxyz ; ii++ ){ 00458 val = CABS(cfim[ii]) ; 00459 if( val > bot && val < top ) cfim[ii].r = cfim[ii].i = 0.0 ; 00460 } 00461 } 00462 break ; 00463 } 00464 } 00465 00466 /*----- apply threshold? -----*/ 00467 00468 if( edit_thresh > 0.0 && iv_thr >= 0 ){ 00469 #ifdef AFNI_DEBUG 00470 int nthresh = 0 ; 00471 # define THADD (nthresh++) 00472 #else 00473 # define THADD /* nada */ 00474 #endif 00475 00476 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: apply threshold\n") ; 00477 00478 switch( thr_type ){ 00479 00480 /** threshold datum is shorts **/ 00481 00482 case MRI_short:{ 00483 short thrplu , thrmin ; 00484 float fplu = edit_thresh / thrfac ; 00485 if( fplu > 32767.0 ){ 00486 fprintf(stderr,"\n*** -1thresh out of range: reset to %g\n", 00487 32767.0 * thrfac ) ; 00488 fplu = 32767.0 ; 00489 } 00490 thrplu = (short) fplu ; 00491 thrmin = -thrplu ; 00492 #ifdef AFNI_DEBUG 00493 { char str[256] ; 00494 sprintf(str,"short threshold = %d\n",(int)thrplu) ; STATUS(str) ; } 00495 #endif 00496 switch( fim_type ){ 00497 case MRI_short: /* fim datum is shorts */ 00498 for( ii=0 ; ii < nxyz ; ii++ ) 00499 if( sthr[ii] < thrplu && sthr[ii] > thrmin ){ sfim[ii] = 0 ; THADD ; } 00500 break ; 00501 00502 case MRI_byte: /* fim datum is bytes */ 00503 for( ii=0 ; ii < nxyz ; ii++ ) 00504 if( sthr[ii] < thrplu && sthr[ii] > thrmin ){ bfim[ii] = 0 ; THADD ; } 00505 break ; 00506 00507 case MRI_float: /* fim datum is floats */ 00508 for( ii=0 ; ii < nxyz ; ii++ ) 00509 if( sthr[ii] < thrplu && sthr[ii] > thrmin ){ ffim[ii] = 0.0 ; THADD ; } 00510 break ; 00511 00512 case MRI_complex: /* fim datum is complex */ 00513 for( ii=0 ; ii < nxyz ; ii++ ) 00514 if( sthr[ii] < thrplu && sthr[ii] > thrmin ){ 00515 cfim[ii].r = cfim[ii].i = 0.0 ; THADD ; 00516 } 00517 break ; 00518 } 00519 } 00520 break ; 00521 00522 /** threshold datum is bytes **/ 00523 00524 case MRI_byte:{ 00525 byte thrplu ; 00526 float fplu = edit_thresh / thrfac ; 00527 if( fplu > 255.0 ){ 00528 fprintf(stderr,"\n*** -1thresh out of range: reset to %g\n", 00529 255.0 * thrfac ) ; 00530 fplu = 255.0 ; 00531 } 00532 thrplu = (byte) fplu ; 00533 #ifdef AFNI_DEBUG 00534 { char str[256] ; 00535 sprintf(str,"byte threshold = %d\n",(int)thrplu) ; STATUS(str) ; } 00536 #endif 00537 switch( fim_type ){ 00538 case MRI_short: /* fim datum is shorts */ 00539 for( ii=0 ; ii < nxyz ; ii++ ) if( bthr[ii] < thrplu ){ sfim[ii] = 0 ; THADD ; } 00540 break ; 00541 00542 case MRI_byte: /* fim datum is bytes */ 00543 for( ii=0 ; ii < nxyz ; ii++ ) if( bthr[ii] < thrplu ){ bfim[ii] = 0 ; THADD ; } 00544 break ; 00545 00546 case MRI_float: /* fim datum is floats */ 00547 for( ii=0 ; ii < nxyz ; ii++ ) if( bthr[ii] < thrplu ){ ffim[ii] = 0.0 ; THADD ; } 00548 break ; 00549 00550 case MRI_complex: /* fim datum is complex */ 00551 for( ii=0 ; ii < nxyz ; ii++ ) 00552 if( bthr[ii] < thrplu ){ 00553 cfim[ii].r = cfim[ii].i = 0.0 ; THADD ; 00554 } 00555 break ; 00556 } 00557 } 00558 break ; 00559 00560 /** threshold datum is floats **/ 00561 00562 case MRI_float:{ 00563 float thrplu , thrmin ; 00564 thrplu = edit_thresh ; if( thrfac > 0.0 ) thrplu /= thrfac ; 00565 thrmin = -thrplu ; 00566 #ifdef AFNI_DEBUG 00567 { char str[256] ; 00568 sprintf(str,"float threshold = %g\n",thrplu) ; STATUS(str) ; } 00569 #endif 00570 switch( fim_type ){ 00571 case MRI_short: /* fim datum is shorts */ 00572 for( ii=0 ; ii < nxyz ; ii++ ) 00573 if( fthr[ii] < thrplu && fthr[ii] > thrmin ){ sfim[ii] = 0 ; THADD ; } 00574 break ; 00575 00576 case MRI_byte: /* fim datum is bytes */ 00577 for( ii=0 ; ii < nxyz ; ii++ ) 00578 if( fthr[ii] < thrplu && fthr[ii] > thrmin ){ bfim[ii] = 0 ; THADD ; } 00579 break ; 00580 00581 case MRI_float: /* fim datum is floats */ 00582 for( ii=0 ; ii < nxyz ; ii++ ) 00583 if( fthr[ii] < thrplu && fthr[ii] > thrmin ){ ffim[ii] = 0.0 ; THADD ; } 00584 break ; 00585 00586 case MRI_complex: /* fim datum is complex */ 00587 for( ii=0 ; ii < nxyz ; ii++ ) 00588 if( fthr[ii] < thrplu && fthr[ii] > thrmin ){ 00589 cfim[ii].r = cfim[ii].i = 0.0 ; THADD ; 00590 } 00591 break ; 00592 } 00593 } 00594 break ; 00595 } 00596 #ifdef AFNI_DEBUG 00597 { char str[256] ; 00598 sprintf(str,"number thresholded to zero = %d",nthresh) ; 00599 STATUS(str) ; } 00600 #endif 00601 } 00602 00603 /*----- blur? -----*/ 00604 00605 if( edit_blur > 0.0 ){ 00606 00607 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: blurring fim\n") ; 00608 00609 EDIT_blur_volume( nx,ny,nz, dx,dy,dz , fim_type,vfim , edit_blur ) ; 00610 } 00611 00612 /*----- threshold blur? -----*/ /* 4 Oct 1996 */ 00613 if(( edit_thrblur > 0.0) && (vthr != NULL) ){ 00614 00615 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: blurring threshold\n") ; 00616 00617 EDIT_blur_volume( nx,ny,nz, dx,dy,dz , thr_type,vthr , edit_thrblur ) ; 00618 } 00619 00620 00621 /*----- zvol? -----*/ 00622 00623 if( edit_zvol ){ 00624 THD_ivec3 iv1 , iv2 ; 00625 int ix1,ix2 , jy1,jy2 , kz1,kz2 , jj,kk ; 00626 00627 iv1 = THD_3dmm_to_3dind(dset,TEMP_FVEC3(edopt->zv_x1,edopt->zv_y1,edopt->zv_z1)); 00628 iv2 = THD_3dmm_to_3dind(dset,TEMP_FVEC3(edopt->zv_x2,edopt->zv_y2,edopt->zv_z2)); 00629 00630 ix1 = iv1.ijk[0] ; ix2 = iv2.ijk[0] ; 00631 jy1 = iv1.ijk[1] ; jy2 = iv2.ijk[1] ; 00632 kz1 = iv1.ijk[2] ; kz2 = iv2.ijk[2] ; 00633 00634 if( ix1 > ix2 ){ ii=ix1 ; ix1=ix2 ; ix2=ii ; } 00635 if( jy1 > jy2 ){ ii=jy1 ; jy1=jy2 ; jy2=ii ; } 00636 if( kz1 > kz2 ){ ii=kz1 ; kz1=kz2 ; kz2=ii ; } 00637 00638 #ifdef AFNI_DEBUG 00639 { char str[256] ; 00640 sprintf(str,"edit_zvol: x1=%g x2=%g y1=%g y2=%g z1=%g z2=%g", 00641 edopt->zv_x1,edopt->zv_x2,edopt->zv_y1,edopt->zv_y2,edopt->zv_z1,edopt->zv_z2) ; 00642 STATUS(str) ; 00643 sprintf(str," : ix1=%d ix2=%d jy1=%d jy2=%d kz1=%d kz2=%d", 00644 ix1,ix2,jy1,jy2,kz1,kz2) ; 00645 STATUS(str) ; } 00646 #endif 00647 00648 if( verbose ) 00649 fprintf(stderr,"--- EDIT_one_dataset: zeroing indexes [%d,%d]x[%d,%d]x[%d,%d]\n", 00650 ix1,ix2,jy1,jy2,kz1,kz2 ) ; 00651 00652 for( kk=kz1 ; kk <= kz2 ; kk++ ){ 00653 for( jj=jy1 ; jj <= jy2 ; jj++ ){ 00654 switch( fim_type ){ 00655 case MRI_short: 00656 for( ii=ix1 ; ii <= ix2 ; ii++ ) sfim[ii+jj*nx+kk*nxy] = 0 ; 00657 break ; 00658 00659 case MRI_byte: 00660 for( ii=ix1 ; ii <= ix2 ; ii++ ) bfim[ii+jj*nx+kk*nxy] = 0 ; 00661 break ; 00662 00663 case MRI_float: 00664 for( ii=ix1 ; ii <= ix2 ; ii++ ) ffim[ii+jj*nx+kk*nxy] = 0 ; 00665 break ; 00666 00667 case MRI_complex: 00668 for( ii=ix1 ; ii <= ix2 ; ii++ ) 00669 cfim[ii+jj*nx+kk*nxy].r = cfim[ii+jj*nx+kk*nxy].i = 0 ; 00670 break ; 00671 } 00672 } 00673 } 00674 } 00675 00676 /*----- form clusters? -----*/ 00677 00678 rmm = clust_rmm ; 00679 vmul = clust_vmul ; 00680 00681 if( rmm >= 0.0 ){ /* do clustering? */ 00682 00683 MCW_cluster_array * clbig ; 00684 MCW_cluster * cl ; 00685 00686 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: clustering with rmm=%g vmul=%g\n", 00687 rmm,vmul ) ; 00688 00689 /*----- Erosion and dilation of clusters -----*/ /* 17 June 1998 */ 00690 if (erode_pv > 0.0) 00691 MCW_erode_clusters (nx, ny, nz, dx, dy, dz, fim_type, vfim, rmm, 00692 erode_pv, dilate); 00693 00694 00695 STATUS("clustering") ; 00696 00697 if( vmul >= 0.0 ) 00698 ptmin = (int)( vmul / dxyz + 0.99 ) ; 00699 else 00700 ptmin = (int) fabs(vmul) ; /* 30 Apr 2002 */ 00701 00702 vmul = MAX(1,ptmin) * dxyz ; /* for use below */ 00703 00704 clar = MCW_find_clusters( nx,ny,nz , dx,dy,dz , fim_type,vfim , rmm ) ; 00705 nclu = 0 ; 00706 00707 if( clar != NULL ){ 00708 INIT_CLARR(clbig) ; 00709 for( iclu=0 ; iclu < clar->num_clu ; iclu++ ){ 00710 cl = clar->clar[iclu] ; 00711 if( cl->num_pt >= ptmin ){ /* big enough */ 00712 ADDTO_CLARR(clbig,cl) ; /* copy pointer */ 00713 clar->clar[iclu] = NULL ; /* null out original */ 00714 nclu++ ; 00715 } 00716 } 00717 DESTROY_CLARR(clar) ; 00718 clar = clbig ; 00719 if( nclu == 0 || clar == NULL || clar->num_clu == 0 ){ 00720 printf("*** NO CLUSTERS FOUND ***\n") ; 00721 if( clar != NULL ) DESTROY_CLARR(clar) ; 00722 EXRETURN ; 00723 } 00724 SORT_CLARR(clar) ; 00725 } 00726 00727 if( nclu == 0 ){ /* no data left */ 00728 STATUS("no data left after cluster edit!") ; 00729 DESTROY_CLARR(clar) ; 00730 EXRETURN ; 00731 } 00732 00733 #ifdef AFNI_DEBUG 00734 { char str[256] ; 00735 sprintf(str,"number clusters = %d",nclu) ; STATUS(str) ; } 00736 #endif 00737 00738 /*----- edit clusters? -----*/ /* 10 Sept 1996 */ 00739 if (edit_clust > ECFLAG_SAME) 00740 EDIT_cluster_array (clar, edit_clust, dxyz, vmul); 00741 if (edit_clust == ECFLAG_SIZE || edit_clust == ECFLAG_ORDER) 00742 DSET_BRICK_FACTOR(dset,iv_fim) = 0.0; 00743 00744 for( iclu=0 ; iclu < clar->num_clu ; iclu++ ) 00745 if( clar->clar[iclu] != NULL && clar->clar[iclu]->num_pt > 0 ){ 00746 MCW_cluster_to_vol( nx,ny,nz , fim_type,vfim , clar->clar[iclu] ) ; 00747 } else { 00748 } 00749 00750 DESTROY_CLARR(clar) ; 00751 } 00752 00753 00754 /*----- filter? -----*/ /* 11 Sept 1996 */ 00755 if (filter_opt > FCFLAG_NONE){ 00756 00757 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: filtering fim\n") ; 00758 00759 EDIT_filter_volume (nx, ny, nz, dx, dy, dz, fim_type, vfim, 00760 filter_opt, filter_rmm , edopt->fmask , edopt->fexpr ); 00761 } 00762 00763 00764 /*----- threshold filter? -----*/ /* 1 Oct 1996 */ 00765 if ((thrfilter_opt > FCFLAG_NONE) && (vthr != NULL)){ 00766 00767 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: filtering thr\n") ; 00768 00769 EDIT_filter_volume (nx, ny, nz, dx, dy, dz, thr_type, vthr, 00770 thrfilter_opt, thrfilter_rmm , edopt->fmask , edopt->fexpr ); 00771 } 00772 00773 00774 /*----- scale? -----*/ 00775 00776 #ifdef ALLOW_SCALE_TO_MAX 00777 if( edit_scale ){ 00778 STATUS("scale") ; 00779 00780 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: scaling fim to max\n") ; 00781 00782 MCW_scale_to_max( nx,ny,nz , fim_type , vfim ) ; 00783 } 00784 #endif 00785 00786 /*----- mult? -----*/ 00787 /*--- correction for scaling of short and byte bricks (13 Sept. 1996) ---*/ 00788 00789 if( edit_mult != 0.0 ){ 00790 STATUS("mult") ; 00791 00792 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: multiplying fim\n") ; 00793 00794 switch( fim_type ){ 00795 case MRI_short: 00796 if (fimfac > 0) 00797 DSET_BRICK_FACTOR(dset,iv_fim) = 00798 DSET_BRICK_FACTOR(dset,iv_fim) * edit_mult ; 00799 else 00800 for( ii=0 ; ii < nxyz ; ii++ ) sfim[ii] *= edit_mult ; 00801 break ; 00802 00803 case MRI_byte : 00804 if (fimfac > 0) 00805 DSET_BRICK_FACTOR(dset,iv_fim) = 00806 DSET_BRICK_FACTOR(dset,iv_fim) * edit_mult ; 00807 else 00808 for( ii=0 ; ii < nxyz ; ii++ ) bfim[ii] *= edit_mult ; 00809 break ; 00810 00811 case MRI_float: for( ii=0 ; ii < nxyz ; ii++ ) ffim[ii] *= edit_mult ; 00812 break ; 00813 00814 case MRI_complex: for( ii=0 ; ii < nxyz ; ii++ ) 00815 cfim[ii].r *= edit_mult , cfim[ii].i *= edit_mult ; 00816 break ; 00817 } 00818 } 00819 00820 /*----- 17 Sep 1998: conversion to z-score? -----*/ 00821 00822 if( edopt->zscore ){ /* 07 Jun 1999! How did this get lost? */ 00823 int kv = DSET_BRICK_STATCODE(dset,iv_fim) ; 00824 float par[2] ; 00825 00826 if( FUNC_IS_STAT(kv) && kv != FUNC_ZT_TYPE ){ 00827 00828 #if 0 00829 fprintf(stderr," -1zscore: converting\n") ; 00830 #endif 00831 00832 if( verbose ) fprintf(stderr,"--- EDIT_one_dataset: converting to zscore\n") ; 00833 00834 EDIT_zscore_vol( nxyz , fim_type , fimfac , vfim , 00835 kv , DSET_BRICK_STATAUX(dset,iv_fim) ) ; 00836 00837 if( ISBUCKET(dset) ){ 00838 00839 #if 0 00840 fprintf(stderr," -1zscore: bucketing\n") ; 00841 #endif 00842 00843 par[0] = FUNC_ZT_TYPE ; 00844 par[1] = 0 ; 00845 EDIT_dset_items( dset , ADN_brick_stataux_one+iv_fim,par , ADN_none ) ; 00846 00847 } else if( ISFUNC(dset) && 00848 FUNC_IS_STAT(dset->func_type) && 00849 iv_fim == FUNC_ival_thr[dset->func_type] ){ 00850 00851 #if 0 00852 fprintf(stderr," -1zscore: retyping\n") ; 00853 #endif 00854 00855 dset->func_type = FUNC_ZT_TYPE ; 00856 dset->stat_aux[0] = 0.0 ; 00857 00858 } else { 00859 fprintf(stderr,"*** -1zscore error: non-bucket & non-func!\n") ; 00860 } 00861 00862 if( fim_type == MRI_short ) 00863 DSET_BRICK_FACTOR(dset,iv_fim) = 1.0 / FUNC_ZT_SCALE_SHORT ; 00864 } 00865 } 00866 00867 /*------ DONE! -----*/ 00868 00869 EXRETURN ; 00870 } |