/**************************************************************************** ***** This file is mean to be #include-d, especially for 3dXClustSim.c ***** ****************************************************************************/ #include "mrilib.h" /* get thread number if compiled with OpenMP -- or zero if not */ #undef DECLARE_ithr /* 30 Nov 2015 */ #ifdef USE_OMP # define DECLARE_ithr const int ithr=omp_get_thread_num() #else # define DECLARE_ithr const int ithr=0 #endif /*---------------------------------------------------------------------------*/ /* Cluster definition. Index type ind_t can be byte or short. */ #define USE_UBYTE /*====================================*/ #ifndef IND_T #ifdef USE_UBYTE /* for grids <= 255 on each side */ # define DALL 128 # define MAX_IND 255u typedef unsigned char ind_t ; # define IND_T unsigned char #else /*==== for grids <= 65535 on each side ====*/ # define DALL 512 # define MAX_IND 32767u typedef unsigned short ind_t ; # define IND_T unsigned short #endif #endif /*====================================*/ /*----- struct to hold a single cluster of points -----*/ typedef struct { int npt , /* number of points assigned */ nall , /* number of points allocated */ norig , /* number of original pts (before dilation) */ nbcount ; /* for dilation decisions (cf. 3dXClustSim) */ float fomh[3] ; /* Figure Of Merit for hpow=0, 1, 2 */ float cth[3] ; /* FOM threshold for hpow=0, 1, 2 */ int hmask ; /* which h FOM threshs were passed */ ind_t *ip , *jp , *kp ; /* 3D indexes for each point */ int *ijk ; /* 1D index for each point */ int ijkmin,ijkmax ; /* smallest and largest ijk[] values */ int indx ; /* some code for 3dXClustSim */ } Xcluster ; #define MIN_CLUST 5 /* smallest cluster size allowed (voxels) */ static int min_clust = MIN_CLUST ; /* 21 Sep 2017 */ void set_Xmin_clust( int mc ){ min_clust = MAX(mc,2); return; } /*----- struct to hold a bunch of clusters -----*/ typedef struct { int nclu , nall ; Xcluster **xclu ; } Xcluster_array ; /*----------------------------------------------------------------------------*/ #define CREATE_Xcluster_array(xcar,siz) \ do{ (xcar) = (Xcluster_array *)malloc(sizeof(Xcluster_array)) ; \ (xcar)->nclu = 0 ; xcar->nall = (siz) ; \ (xcar)->xclu = (Xcluster **)calloc(sizeof(Xcluster *),(siz)) ; \ } while(0) #define ADDTO_Xcluster_array(xcar,xc) \ do{ if( (xcar)->nclu >= (xcar)->nall ){ \ (xcar)->nall = (xcar)->nclu + 32 ; \ (xcar)->xclu = (Xcluster **)realloc((xcar)->xclu, \ sizeof(Xcluster **)*(xcar)->nall) ; \ } \ (xcar)->xclu[(xcar)->nclu++] = (xc) ; \ } while(0) #define DESTROY_Xcluster_array(xcar) \ do{ int cc ; \ for( cc=0 ; cc < (xcar)->nclu ; cc++ ) \ DESTROY_Xcluster((xcar)->xclu[cc]) ; \ free((xcar)->xclu) ; free(xcar) ; \ } while(0) #define MERGE_Xcluster_arrays(xcar,ycar) \ do{ if( (xcar)->nclu+(ycar)->nclu > (xcar)->nall ){ \ (xcar)->nall = (xcar)->nclu+(ycar)->nclu ; \ (xcar)->xclu = (Xcluster **)realloc((xcar)->xclu, \ sizeof(Xcluster **)*(xcar)->nall) ; \ } \ memcpy( (xcar)->xclu+(xcar)->nclu , \ (ycar)->xclu , sizeof(Xcluster *)*(ycar)->nclu ) ; \ free((ycar)->xclu) ; free(ycar) ; \ } while(0) /*----------------------------------------------------------------------------*/ /* Create a cluster with initial array allocation of siz */ #define CREATE_Xcluster(xc,siz) \ do{ xc = (Xcluster *)malloc(sizeof(Xcluster)) ; \ xc->npt = xc->norig = xc->nbcount = 0 ; xc->nall = (siz) ; \ xc->fomh[0] = xc->fomh[1] = xc->fomh[2] = 0.0f ; \ xc->cth [0] = xc->cth [1] = xc->cth [2] = 0.0f ; \ xc->ip = (ind_t *)malloc(sizeof(ind_t)*(siz)) ; \ xc->jp = (ind_t *)malloc(sizeof(ind_t)*(siz)) ; \ xc->kp = (ind_t *)malloc(sizeof(ind_t)*(siz)) ; \ xc->ijk = (int *) malloc(sizeof(int) *(siz)) ; \ xc->ijkmin = xc->ijkmax = xc->indx = -666 ; \ } while(0) #define DESTROY_Xcluster(xc) \ do{ if( (xc) != NULL ){ \ free((xc)->ip); free((xc)->jp); free((xc)->kp); \ free((xc)->ijk); free(xc); \ }} while(0) #if 0 /************************************************************************/ /*----------------------------------------------------------------------------*/ /* Copy one cluster's data over another's */ void copyover_Xcluster( Xcluster *xcin , Xcluster *xcout ) { int nin ; if( xcin == NULL || xcout == NULL ) return ; nin = xcin->npt ; if( nin > xcout->nall ){ xcout->nall = nin ; xcout->ip = (ind_t *)realloc(xcout->ip ,sizeof(ind_t)*nin) ; xcout->jp = (ind_t *)realloc(xcout->jp ,sizeof(ind_t)*nin) ; xcout->kp = (ind_t *)realloc(xcout->kp ,sizeof(ind_t)*nin) ; xcout->ijk = (int *) realloc(xcout->ijk,sizeof(int) *nin) ; } AA_memcpy( xcout->ip , xcin->ip , sizeof(ind_t)*nin ) ; AA_memcpy( xcout->jp , xcin->jp , sizeof(ind_t)*nin ) ; AA_memcpy( xcout->kp , xcin->kp , sizeof(ind_t)*nin ) ; AA_memcpy( xcout->ijk, xcin->ijk, sizeof(int) *nin ) ; xcout->fomh[0] = xcin->fomh[0] ; xcout->fomh[1] = xcin->fomh[1] ; xcout->fomh[2] = xcin->fomh[2] ; xcout->cth [0] = xcin->cth [0] ; xcout->cth [1] = xcin->cth [1] ; xcout->cth [2] = xcin->cth [2] ; xcout->npt = nin ; xcout->norig = xcin->norig ; xcout->nbcount = xcin->nbcount ; xcout->ijkmin = xcin->ijkmin ; xcout->ijkmax = xcin->ijkmax ; xcout->indx = xcin->indx ; return ; } /*----------------------------------------------------------------------------*/ /* Create a new cluster that is a copy of the input */ Xcluster * copy_Xcluster( Xcluster *xcc ) { Xcluster *xccout ; if( xcc == NULL ) return NULL ; CREATE_Xcluster(xccout,xcc->npt) ; copyover_Xcluster( xcc , xccout ) ; return xccout ; } #endif /***********************************************************************/ /*----------------------------------------------------------------------------*/ /* Add a point to current cluster (if far is nonzero at this point). For use only in the function directly below! *//*--------------------------------------------------------------------------*/ static float **cthar0 = NULL ; /* arrays to keep track of cluster */ static int *ncthar0 = NULL ; /* thresholds at each point and in */ static int *kcthar0 = NULL ; /* each OpenMP thread -- */ /* see mri_multi_threshold_setup() */ static float **cthar1 = NULL ; static int *ncthar1 = NULL ; static int *kcthar1 = NULL ; static float **cthar2 = NULL ; static int *ncthar2 = NULL ; static int *kcthar2 = NULL ; static int cth_mode = 2 ; /* 0 = mean , 1 = median, 2 = cth_perc% */ static float cth_perc = 95.0f ; #define ADDTO_CTHAR0(val,ith) \ do{ if( kcthar0[ith] >= ncthar0[ith] ){ \ ncthar0[ith] = 2*kcthar0[ith] ; \ cthar0[ith] = (float *)realloc(cthar0[ith],sizeof(float)*ncthar0[ith]); \ } \ cthar0[ith][kcthar0[ith]++] = (val) ; \ } while(0) #define ADDTO_CTHAR1(val,ith) \ do{ if( kcthar1[ith] >= ncthar1[ith] ){ \ ncthar1[ith] = 2*kcthar1[ith] ; \ cthar1[ith] = (float *)realloc(cthar1[ith],sizeof(float)*ncthar1[ith]); \ } \ cthar1[ith][kcthar1[ith]++] = (val) ; \ } while(0) #define ADDTO_CTHAR2(val,ith) \ do{ if( kcthar2[ith] >= ncthar2[ith] ){ \ ncthar2[ith] = 2*kcthar2[ith] ; \ cthar2[ith] = (float *)realloc(cthar2[ith],sizeof(float)*ncthar2[ith]); \ } \ cthar2[ith][kcthar2[ith]++] = (val) ; \ } while(0) #define CPUT_point(i,j,k) \ do{ int pqr = (i)+(j)*nx+(k)*nxy , npt=(xcc)->npt ; \ if( far[pqr] != 0.0f ){ \ if( npt >= (xcc)->nall ){ \ (xcc)->nall = npt + DALL + (xcc)->nall/2 ; \ (xcc)->ip = (ind_t *)realloc((xcc)->ip,sizeof(ind_t)*(xcc)->nall) ; \ (xcc)->jp = (ind_t *)realloc((xcc)->jp,sizeof(ind_t)*(xcc)->nall) ; \ (xcc)->kp = (ind_t *)realloc((xcc)->kp,sizeof(ind_t)*(xcc)->nall) ; \ (xcc)->ijk= (int *) realloc((xcc)->ijk,sizeof(int) *(xcc)->nall) ; \ } \ (xcc)->ip[npt] = (i); (xcc)->jp[npt] = (j); (xcc)->kp[npt] = (k); \ (xcc)->ijk[npt] = pqr ; \ if( (xcc)->ijkmin > pqr ) (xcc)->ijkmin = pqr ; \ if( (xcc)->ijkmax < pqr ) (xcc)->ijkmax = pqr ; \ (xcc)->npt++ ; (xcc)->norig++ ; \ (xcc)->fomh[0] += 1.0f ; \ (xcc)->fomh[1] += fabsf(far[pqr]) ; \ (xcc)->fomh[2] += far[pqr]*far[pqr] ; \ if( docim0 ) ADDTO_CTHAR0(car0[pqr],ithr) ; \ if( docim1 ) ADDTO_CTHAR1(car1[pqr],ithr) ; \ if( docim2 ) ADDTO_CTHAR2(car2[pqr],ithr) ; \ far[pqr] = 0.0f ; \ } } while(0) /*----------------------------------------------------------------------------*/ /* Find clusters of nonzero voxels: fim = thresholded float image (will be zero-ed out at end) nnlev = 1 or 2 or 3 (clustering neighborliness) thtype = threshold type (0 = float *, 1 = MRI_IMAGE *) thptrX = pointer to threshold (a single float, or an image) *//*--------------------------------------------------------------------------*/ Xcluster_array * find_Xcluster_array( MRI_IMAGE *fim, int nnlev, int indx , int thtype , void *thptr0 , void *thptr1 , void *thptr2 ) { Xcluster *xcc=NULL ; Xcluster_array *xcar=NULL ; float *far,*car0=NULL,*car1=NULL,*car2=NULL ; int ii,jj,kk, icl , ijk , ijk_last ; int ip,jp,kp , im,jm,km , nx,ny,nz,nxy,nxyz ; const int do_nn2=(nnlev > 1) , do_nn3=(nnlev > 2) ; int docim0=0, docim1=0, docim2=0, doccc=0 ; ; MRI_IMAGE *cim0=NULL, *cim1=NULL, *cim2=NULL ; float *fpt0=NULL, *fpt1=NULL, *fpt2=NULL ; float fth0=0.0f, fth1=0.0f, fth2=0.0f ; int dofth0=0 ,dofth1=0 ,dofth2=0 , dofff=0 ; int dott0 =0 ,dott1 =0 ,dott2 =0 , dottt=0 ; DECLARE_ithr ; far = MRI_FLOAT_PTR(fim) ; nx = fim->nx; ny = fim->ny; nxy = nx*ny; nz = fim->nz; nxyz = nxy*nz; /* setup thresholds: spatially constant (thtype==0) or variable (==1) */ switch( thtype ){ case 0: fpt0 = (float *)thptr0; if( fpt0 != NULL ){ fth0 = *fpt0; dofth0 = (fth0 > 0.0f); } fpt1 = (float *)thptr1; if( fpt1 != NULL ){ fth1 = *fpt1; dofth1 = (fth1 > 0.0f); } fpt2 = (float *)thptr2; if( fpt2 != NULL ){ fth2 = *fpt2; dofth2 = (fth2 > 0.0f); } /* INFO_message("find_Xcluster_array: fth0=%.1f fth1=%.1f fth2=%.1f",fth0,fth1,fth2) ; */ dofff = (dofth0||dofth1||dofth2) ; break ; case 1: cim0 = (MRI_IMAGE *)thptr0; docim0 = (cim0!=NULL); if( docim0 ) car0 = MRI_FLOAT_PTR(cim0); cim1 = (MRI_IMAGE *)thptr1; docim1 = (cim1!=NULL); if( docim1 ) car1 = MRI_FLOAT_PTR(cim1); cim2 = (MRI_IMAGE *)thptr2; docim2 = (cim2!=NULL); if( docim2 ) car2 = MRI_FLOAT_PTR(cim2); doccc =(docim0||docim1||docim2) ; break ; default: doccc = dofff = 0 ; if( thptr0 != NULL || thptr1 != NULL || thptr2 != NULL ) WARNING_message("bad thtype=%d in find_Xcluster_array :(",thtype) ; break ; } dott0 = (dofth0 || docim0) ; dott1 = (dofth1 || docim1) ; dott2 = (dofth2 || docim2) ; dottt = (doccc || dofff) ; ijk_last = 0 ; /* start scanning at the {..wait for it..} start */ while(1){ /* find next nonzero point in far array (not far away) */ for( ijk=ijk_last ; ijk < nxyz ; ijk++ ) if( far[ijk] != 0.0f ) break ; if( ijk == nxyz ) break ; /* didn't find any! */ ijk_last = ijk+1 ; /* start here next time */ IJK_TO_THREE(ijk, ii,jj,kk , nx,nxy) ; /* 3D coords of this point */ /* build a new cluster starting with this 1 point */ if( xcc == NULL ) CREATE_Xcluster(xcc,16) ; /* initialize to have just 16 potential points */ xcc->ip[0] = (ind_t)ii; xcc->jp[0] = (ind_t)jj; xcc->kp[0] = (ind_t)kk; xcc->ijk[0] = xcc->ijkmin = xcc->ijkmax = ijk; xcc->indx = indx ; xcc->npt = xcc->norig = 1 ; xcc->hmask = 0 ; xcc->fomh[0] = 1.0f ; xcc->fomh[1] = fabsf(far[ijk]) ; xcc->fomh[2] = far[ijk]*far[ijk] ; far[ijk] = 0.0f ; xcc->cth[0] = xcc->cth[1] = xcc->cth[2] = 0.0f ; if( docim0 ){ cthar0[ithr][0] = car0[ijk]; kcthar0[ithr] = 1; } else { kcthar0[ithr] = 0; } if( docim1 ){ cthar1[ithr][0] = car1[ijk]; kcthar1[ithr] = 1; } else { kcthar1[ithr] = 0; } if( docim2 ){ cthar2[ithr][0] = car2[ijk]; kcthar2[ithr] = 1; } else { kcthar2[ithr] = 0; } /* loop over points in cluster, checking their neighbors, growing the cluster if we find any that belong therein */ for( icl=0 ; icl < xcc->npt ; icl++ ){ ii = xcc->ip[icl]; jj = xcc->jp[icl]; kk = xcc->kp[icl]; im = ii-1 ; jm = jj-1 ; km = kk-1 ; /* minus 1 indexes */ ip = ii+1 ; jp = jj+1 ; kp = kk+1 ; /* plus 1 indexes */ /* CPUT_point(i,j,k) only does something if far[i,j,k] is nonzero */ if( im >= 0 ){ CPUT_point(im,jj,kk) ; /* 1NN */ if( do_nn2 ){ if( jm >= 0 ) CPUT_point(im,jm,kk) ; /* 2NN */ if( jp < ny ) CPUT_point(im,jp,kk) ; /* 2NN */ if( km >= 0 ) CPUT_point(im,jj,km) ; /* 2NN */ if( kp < nz ) CPUT_point(im,jj,kp) ; /* 2NN */ if( do_nn3 ){ if( jm >= 0 && km >= 0 ) CPUT_point(im,jm,km) ; /* 3NN */ if( jm >= 0 && kp < nz ) CPUT_point(im,jm,kp) ; /* 3NN */ if( jp < ny && km >= 0 ) CPUT_point(im,jp,km) ; /* 3NN */ if( jp < ny && kp < nz ) CPUT_point(im,jp,kp) ; /* 3NN */ }}} if( ip < nx ){ CPUT_point(ip,jj,kk) ; /* 1NN */ if( do_nn2 ){ if( jm >= 0 ) CPUT_point(ip,jm,kk) ; /* 2NN */ if( jp < ny ) CPUT_point(ip,jp,kk) ; /* 2NN */ if( km >= 0 ) CPUT_point(ip,jj,km) ; /* 2NN */ if( kp < nz ) CPUT_point(ip,jj,kp) ; /* 2NN */ if( do_nn3 ){ if( jm >= 0 && km >= 0 ) CPUT_point(ip,jm,km) ; /* 3NN */ if( jm >= 0 && kp < nz ) CPUT_point(ip,jm,kp) ; /* 3NN */ if( jp < ny && km >= 0 ) CPUT_point(ip,jp,km) ; /* 3NN */ if( jp < ny && kp < nz ) CPUT_point(ip,jp,kp) ; /* 3NN */ }}} if( jm >= 0 ){ CPUT_point(ii,jm,kk) ; /* 1NN */ if( do_nn2 ){ if( km >= 0 ) CPUT_point(ii,jm,km) ; /* 2NN */ if( kp < nz ) CPUT_point(ii,jm,kp) ; /* 2NN */ }} if( jp < ny ){ CPUT_point(ii,jp,kk) ; /* 1NN */ if( do_nn2 ){ if( km >= 0 ) CPUT_point(ii,jp,km) ; /* 2NN */ if( kp < nz ) CPUT_point(ii,jp,kp) ; /* 2NN */ }} if( km >= 0 ) CPUT_point(ii,jj,km) ; /* 1NN */ if( kp < nz ) CPUT_point(ii,jj,kp) ; /* 1NN */ } /* since xcc->npt increases when CPUT_point adds a point, the loop continues until finally no new neighbors get added */ if( xcc->npt < min_clust ){ /* too small ==> recycle */ xcc->npt = xcc->norig = 0; continue ; } /* compute FOM thresholds, if they will be used */ if( docim0 ){ /* h=0 */ float cth , qmean=qmean_float(kcthar0[ithr],cthar0[ithr]) , qmed ; switch( cth_mode ){ case 0: cth = qmean; break; case 1: qmed = qmed_float(kcthar0[ithr],cthar0[ithr]); cth = MAX(qmean,qmed); break; default: case 2: qmed = qfrac_float(kcthar0[ithr],0.01f*cth_perc,cthar0[ithr]); cth = MAX(qmean,qmed); break; } xcc->cth[0] = cth ; } else if( dofth0 ) { xcc->cth[0] = fth0 ; } if( docim1 ){ /* h=1 */ float cth , qmean=qmean_float(kcthar1[ithr],cthar1[ithr]) , qmed ; switch( cth_mode ){ case 0: cth = qmean; break; case 1: qmed = qmed_float(kcthar1[ithr],cthar1[ithr]); cth = MAX(qmean,qmed); break; default: case 2: qmed = qfrac_float(kcthar1[ithr],0.01f*cth_perc,cthar1[ithr]); cth = MAX(qmean,qmed); break; } xcc->cth[1] = cth ; } else if( dofth1 ) { xcc->cth[1] = fth1 ; } if( docim2 ){ /* h=2 */ float cth , qmean=qmean_float(kcthar2[ithr],cthar2[ithr]) , qmed ; switch( cth_mode ){ case 0: cth = qmean; break; case 1: qmed = qmed_float(kcthar2[ithr],cthar2[ithr]); cth = MAX(qmean,qmed); break; default: case 2: qmed = qfrac_float(kcthar2[ithr],0.01f*cth_perc,cthar2[ithr]); cth = MAX(qmean,qmed); break; } xcc->cth[2] = cth ; } else if( dofth2 ) { xcc->cth[2] = fth2 ; } /* check thresholds vs cluster FOMs to see if cluster is worthy? */ if( dottt ){ int ngood=0 ; if( dott0 && xcc->fomh[0] >= xcc->cth[0] ) ngood += 1 ; if( dott1 && xcc->fomh[1] >= xcc->cth[1] ) ngood += 2 ; if( dott2 && xcc->fomh[2] >= xcc->cth[2] ) ngood += 4 ; if( ngood == 0 ){ /* didn't pass any threshold ==> recycle */ xcc->npt = xcc->norig = 0; continue; /* cluster is set to be empty */ } else { xcc->hmask = ngood ; } } /* if cluster makes it to here, throw it on the pile */ if( xcar == NULL ) CREATE_Xcluster_array(xcar,8) ; /* create the list */ ADDTO_Xcluster_array(xcar,xcc) ; xcc = NULL ; /* no recycling ==> will create new xcc when needed */ } /* loop until all nonzero points in far[] have been used up */ if( xcc != NULL ) DESTROY_Xcluster(xcc) ; /* exited with unsaved cluster? */ return xcar ; /* could be NULL (if nothing survived the tests) */ } #undef CPUT_point /* not to be used again */ /*----------------------------------------------------------------------------*/ void set_Xcluster_ijkminmax( Xcluster *xcc ) { int ib,it , ii,npt ; if( xcc == NULL ) return ; npt = xcc->npt ; if( npt < 1 ) return ; ib = it = xcc->ijk[0] ; for( ii=1 ; ii < npt ; ii++ ){ if( xcc->ijk[ii] < ib ) ib = xcc->ijk[ii] ; else if( xcc->ijk[ii] > it ) it = xcc->ijk[ii] ; } xcc->ijkmin = ib ; xcc->ijkmax = it ; return ; } /*----------------------------------------------------------------------------*/ /* Allocate the initial cthar arrays for cluster FOM thresholding */ void mri_multi_threshold_setup(void) { int nthr=1 , ithr ; char *eee ; #ifdef USE_OMP nthr = omp_get_max_threads() ; #endif cthar0 = (float **)malloc(sizeof(float *)*nthr) ; ncthar0 = (int * )malloc(sizeof(int)*nthr) ; kcthar0 = (int * )malloc(sizeof(int)*nthr) ; cthar1 = (float **)malloc(sizeof(float *)*nthr) ; ncthar1 = (int * )malloc(sizeof(int)*nthr) ; kcthar1 = (int * )malloc(sizeof(int)*nthr) ; cthar2 = (float **)malloc(sizeof(float *)*nthr) ; ncthar2 = (int * )malloc(sizeof(int)*nthr) ; kcthar2 = (int * )malloc(sizeof(int)*nthr) ; for( ithr=0 ; ithr < nthr ; ithr++ ){ cthar0[ithr] = (float *)malloc(sizeof(float)*4096) ; ncthar0[ithr] = 4096 ; kcthar0[ithr] = 0 ; cthar1[ithr] = (float *)malloc(sizeof(float)*4096) ; ncthar1[ithr] = 4096 ; kcthar1[ithr] = 0 ; cthar2[ithr] = (float *)malloc(sizeof(float)*4096) ; ncthar2[ithr] = 4096 ; kcthar2[ithr] = 0 ; } eee = getenv("AFNI_MTHRESH_MODE") ; if( eee != NULL ){ char *ppp ; cth_mode = 2 ; ppp = strcasestr(eee,"mean") ; if( ppp != NULL ) cth_mode = 0 ; ppp = strcasestr(eee,"median") ; if( ppp != NULL ) cth_mode = 1 ; ppp = strstr(eee,"%") ; if( ppp != NULL && isdigit(*(ppp+1)) ){ cth_perc = (float)strtod(ppp+1,NULL) ; if( cth_perc >= 0.0f && cth_perc <= 100.0f ) cth_mode = 2 ; else { WARNING_message( "AFNI_MTHRESH_MODE says %.1f%% -- out of range 0..100, so using median", cth_perc) ; cth_mode = 1 ; } } } #if 0 switch( cth_mode ){ default: INFO_message("MultiThresh cluster FOM local threshold method = MEAN"); break; case 1: INFO_message("MultiThresh cluster FOM local threshold method = MEDIAN"); break; case 2: INFO_message("MultiThresh cluster FOM local threshold method = cdf %.1f%%",cth_perc); break; } #endif return ; } /*----------------------------------------------------------------------------*/ /* get rid of the cthar arrays */ void mri_multi_threshold_unsetup(void) { int nthr=1 , ithr ; #ifdef USE_OMP nthr = omp_get_max_threads() ; #endif if( cthar0 != NULL ){ for( ithr=0 ; ithr < nthr ; ithr++ ) free(cthar0[ithr]) ; free(cthar0) ; free(ncthar0) ; free(kcthar0) ; } if( cthar1 != NULL ){ for( ithr=0 ; ithr < nthr ; ithr++ ) free(cthar1[ithr]) ; free(cthar1) ; free(ncthar1) ; free(kcthar1) ; } if( cthar2 != NULL ){ for( ithr=0 ; ithr < nthr ; ithr++ ) free(cthar2[ithr]) ; free(cthar2) ; free(ncthar2) ; free(kcthar2) ; } cthar0 = cthar1 = cthar2 = NULL ; ncthar0 = ncthar1 = ncthar2 = NULL ; kcthar0 = kcthar1 = kcthar2 = NULL ; return ; } /*----------------------------------------------------------------------------*/ #define XTHRESH_OUTPUT_MASK 1 #define CIM0(i) ( (cimar0==NULL) ? NULL : IMARR_SUBIM(cimar0,(i)) ) #define CIM1(i) ( (cimar1==NULL) ? NULL : IMARR_SUBIM(cimar1,(i)) ) #define CIM2(i) ( (cimar2==NULL) ? NULL : IMARR_SUBIM(cimar2,(i)) ) /*----------------------------------------------------------------------------*/ /* fim = image to threshold nthr = num thresholds (at least 1) thar = threshold array (length nthr) sid = sideness of threshold (1 or 2) nnlev = NN cluster type (1 or 2 or 3) cimar = array of cluster fom threshold images (if NULL, all clusters >= min_clust voxels are kept) flags = bitwise OR (|) of some or all of the following XTHRESH_OUTPUT_MASK return value will be NULL if nothing was found if nhits!=NULL, the number of 'hits' (nonzero voxels in output) goes there *//*--------------------------------------------------------------------------*/ MRI_IMAGE * mri_multi_threshold_Xcluster_cimar( MRI_IMAGE *fim , int nthr , float *thar , int sid , int nnlev , MRI_IMARR *cimar0 , MRI_IMARR *cimar1 , MRI_IMARR *cimar2 , int flags , int *nhits , MRI_IMARR **allmask ) { MRI_IMAGE *tfim , *qfim=NULL , *cim0,*cim1,*cim2 ; float *tfar , *far , *qfar=NULL , cth,cval,thr ; int ii,nvox , kth , jhp , hm ; Xcluster_array *xcar ; Xcluster *xcc ; int icl,npt=0, *ijkar ; int do_mask = (flags & XTHRESH_OUTPUT_MASK) ; byte *qbyt=NULL ; MRI_IMARR *amask=NULL; MRI_IMAGE *afim=NULL; byte *abyt=NULL; /* bad inputs? */ if( nhits != NULL ) *nhits = 0 ; if( allmask != NULL ) *allmask = NULL ; if( fim == NULL || fim->kind != MRI_float ) return NULL ; if( nthr <= 0 || thar == NULL ) return NULL ; if( cimar0 != NULL && IMARR_COUNT(cimar0) != nthr ) return NULL ; if( cimar1 != NULL && IMARR_COUNT(cimar1) != nthr ) return NULL ; if( cimar2 != NULL && IMARR_COUNT(cimar2) != nthr ) return NULL ; if( allmask != NULL ) INIT_IMARR(amask) ; nvox = fim->nvox ; far = MRI_FLOAT_PTR(fim) ; for( kth=0 ; kth < nthr ; kth++ ){ /* loop over thresholds */ thr = thar[kth] ; cim0 = CIM0(kth) ; cim1 = CIM1(kth) ; cim2 = CIM2(kth) ; tfim = mri_copy(fim) ; tfar = MRI_FLOAT_PTR(tfim) ; /* voxel-wise thresholding */ if( thr != 0.0f ){ if( sid == 2 ){ /*-- 2 sided --*/ thr = fabsf(thr) ; for( ii=0 ; ii < nvox ; ii++ ) if( tfar[ii] > -thr && tfar[ii] < thr ) tfar[ii] = 0.0f ; } else if( thr > 0.0f ){ /*-- 1 sided (positive) --*/ for( ii=0 ; ii < nvox ; ii++ ) if( tfar[ii] < thr ) tfar[ii] = 0.0f ; } else if( thr < 0.0f ){ /*-- 1 sided (negative) --*/ for( ii=0 ; ii < nvox ; ii++ ) if( tfar[ii] > thr ) tfar[ii] = 0.0f ; } } /* clusterize and keep "good" clusters (relative to cim) */ xcar = find_Xcluster_array( tfim , nnlev , 0 , 1,cim0,cim1,cim2 ) ; mri_free(tfim) ; if( xcar == NULL ){ /* we got nuthin at this threshold */ if( amask != NULL ) ADDTO_IMARR(amask,NULL) ; continue ; } /* put "good" clusters into qfim (copying from original input image) */ if( qfim == NULL ){ /* create output image if needed */ if( do_mask ){ qfim = mri_new_conforming(fim,MRI_byte) ; /* zero filled */ qbyt = MRI_BYTE_PTR(qfim) ; } else { qfim = mri_new_conforming(fim,MRI_float) ; /* zero filled */ qfar = MRI_FLOAT_PTR(qfim) ; } } if( amask != NULL ){ afim = mri_new_conforming(fim,MRI_byte) ; /* zero filled */ abyt = MRI_BYTE_PTR(afim) ; } for( icl=0 ; icl < xcar->nclu ; icl++ ){ xcc = xcar->xclu[icl]; npt = xcc->npt; ijkar = xcc->ijk; hm = xcc->hmask; if( do_mask ){ for( ii=0 ; ii < npt ; ii++ ) qbyt[ijkar[ii]]++ ; } else { for( ii=0 ; ii < npt ; ii++ ) qfar[ijkar[ii]] = far[ijkar[ii]] ; } if( amask != NULL ){ for( ii=0 ; ii < npt ; ii++ ) abyt[ijkar[ii]] = hm ; } } DESTROY_Xcluster_array(xcar) ; if( amask != NULL ) ADDTO_IMARR(amask,afim) ; } /* count number of hits */ if( nhits != NULL && qfim != NULL ){ if( do_mask ){ for( npt=ii=0 ; ii < nvox ; ii++ ) npt += (qbyt[ii] != 0) ; } else if( qfar != NULL ){ for( npt=ii=0 ; ii < nvox ; ii++ ) npt += (qfar[ii] != 0.0f) ; } *nhits = npt ; } if( allmask != NULL ) *allmask = amask ; return qfim ; /* nuthin nowhere nohow ==> NULL */ } /*----------------------------------------------------------------------------*/ /* Version of the above with a single per-voxel threshold value and a single cluster-FOM threshold image *//*--------------------------------------------------------------------------*/ MRI_IMAGE * mri_threshold_Xcluster( MRI_IMAGE *fim , float thr , int sid , int nnlev , MRI_IMAGE *cim0 , MRI_IMAGE *cim1 , MRI_IMAGE *cim2 , int flags , int *nhits ) { float thar[1] ; MRI_IMARR *cimar0=NULL,*cimar1=NULL,*cimar2=NULL ; MRI_IMAGE *qfim ; if( cim0 != NULL ){ INIT_IMARR(cimar0); ADDTO_IMARR(cimar0,cim0); } if( cim1 != NULL ){ INIT_IMARR(cimar1); ADDTO_IMARR(cimar1,cim1); } if( cim2 != NULL ){ INIT_IMARR(cimar2); ADDTO_IMARR(cimar2,cim2); } thar[0] = thr ; qfim = mri_multi_threshold_Xcluster_cimar(fim,1,thar,sid,nnlev, cimar0,cimar1,cimar2,flags,nhits,NULL) ; if( cimar0 != NULL ) FREE_IMARR(cimar0) ; if( cimar1 != NULL ) FREE_IMARR(cimar1) ; if( cimar2 != NULL ) FREE_IMARR(cimar2) ; return qfim ; } /*----------------------------------------------------------------------------*/ #define FOM0(i) ( (fomth0==NULL) ? 0.0f : fomth0[(i)] ) #define FOM1(i) ( (fomth1==NULL) ? 0.0f : fomth1[(i)] ) #define FOM2(i) ( (fomth2==NULL) ? 0.0f : fomth2[(i)] ) /*----------------------------------------------------------------------------*/ /* fim = image to threshold nthr = num thresholds (at least 1) thar = threshold array (length nthr) sid = sideness of threshold (1 or 2) nnlev = NN cluster type (1 or 2 or 3) fomth = array of fom threshold values (if all NULL, all clusters >= min_clust voxels are kept) flags = bitwise OR (|) of some or all of the following XTHRESH_OUTPUT_MASK return value will be NULL if nothing was found if nhits!=NULL, the number of 'hits' (nonzero voxels in output) goes there *//*--------------------------------------------------------------------------*/ MRI_IMAGE * mri_multi_threshold_Xcluster_fomth( MRI_IMAGE *fim , int nthr , float *thar , int sid , int nnlev , float *fomth0, float *fomth1, float *fomth2, int flags , int *nhits , MRI_IMARR **allmask ) { MRI_IMAGE *tfim , *qfim=NULL ; float *tfar , *far , *qfar=NULL , cth,cval,thr ; int ii,nvox , kth , jhp , hm ; Xcluster_array *xcar ; Xcluster *xcc ; int icl,npt=0, *ijkar ; int do_mask = (flags & XTHRESH_OUTPUT_MASK) ; byte *qbyt=NULL ; MRI_IMARR *amask=NULL; MRI_IMAGE *afim=NULL; byte *abyt=NULL; float fth0,fth1,fth2 ; /* bad inputs? */ if( nhits != NULL ) *nhits = 0 ; if( allmask != NULL ) *allmask = NULL ; if( fim == NULL || fim->kind != MRI_float ) return NULL ; if( nthr <= 0 || thar == NULL ) return NULL ; if( allmask != NULL ) INIT_IMARR(amask) ; nvox = fim->nvox ; far = MRI_FLOAT_PTR(fim) ; for( kth=0 ; kth < nthr ; kth++ ){ /* loop over thresholds */ thr = thar[kth] ; fth0 = FOM0(kth) ; fth1 = FOM1(kth) ; fth2 = FOM2(kth) ; tfim = mri_copy(fim) ; tfar = MRI_FLOAT_PTR(tfim) ; /* voxel-wise thresholding */ if( thr != 0.0f ){ if( sid == 2 ){ /*-- 2 sided --*/ thr = fabsf(thr) ; for( ii=0 ; ii < nvox ; ii++ ) if( tfar[ii] > -thr && tfar[ii] < thr ) tfar[ii] = 0.0f ; } else if( thr > 0.0f ){ /*-- 1 sided (positive) --*/ for( ii=0 ; ii < nvox ; ii++ ) if( tfar[ii] < thr ) tfar[ii] = 0.0f ; } else if( thr < 0.0f ){ /*-- 1 sided (negative) --*/ for( ii=0 ; ii < nvox ; ii++ ) if( tfar[ii] > thr ) tfar[ii] = 0.0f ; } } /* clusterize and keep "good" clusters (relative to fth) */ xcar = find_Xcluster_array( tfim , nnlev , 0 , 0,&fth0,&fth1,&fth2 ) ; mri_free(tfim) ; if( xcar == NULL ){ /* we got nuthin at this threshold */ if( amask != NULL ) ADDTO_IMARR(amask,NULL) ; continue ; } /* put "good" clusters into qfim (copying from original input image) */ if( qfim == NULL ){ /* create output image if needed */ if( do_mask ){ qfim = mri_new_conforming(fim,MRI_byte) ; /* zero filled */ qbyt = MRI_BYTE_PTR(qfim) ; } else { qfim = mri_new_conforming(fim,MRI_float) ; /* zero filled */ qfar = MRI_FLOAT_PTR(qfim) ; } } if( amask != NULL ){ afim = mri_new_conforming(fim,MRI_byte) ; /* zero filled */ abyt = MRI_BYTE_PTR(afim) ; } for( icl=0 ; icl < xcar->nclu ; icl++ ){ xcc = xcar->xclu[icl]; npt = xcc->npt; ijkar = xcc->ijk; hm = xcc->hmask; if( do_mask ){ for( ii=0 ; ii < npt ; ii++ ) qbyt[ijkar[ii]]++ ; } else { for( ii=0 ; ii < npt ; ii++ ) qfar[ijkar[ii]] = far[ijkar[ii]] ; } if( amask != NULL ){ for( ii=0 ; ii < npt ; ii++ ) abyt[ijkar[ii]] = hm ; } } DESTROY_Xcluster_array(xcar) ; if( amask != NULL ) ADDTO_IMARR(amask,afim) ; } /* count number of hits */ if( nhits != NULL && qfim != NULL ){ if( do_mask ){ for( npt=ii=0 ; ii < nvox ; ii++ ) npt += (qbyt[ii] != 0) ; } else if( qfar != NULL ){ for( npt=ii=0 ; ii < nvox ; ii++ ) npt += (qfar[ii] != 0.0f) ; } *nhits = npt ; } if( allmask != NULL ) *allmask = amask ; return qfim ; /* nuthin nowhere nohow ==> NULL */ }