00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008 #include "parser.h"
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038 void EDIT_filter_volume (int nx, int ny, int nz, float dx, float dy, float dz,
00039 int fim_type, void * vfim, int filter_opt, float filter_rmm,
00040 byte * fmask , char * fexpr )
00041 {
00042 MCW_cluster * mask;
00043 int nxy, nxyz;
00044 int mnum;
00045 int i, j, k, ii, jj, kk;
00046 int ijkvox, ijkma, jma;
00047 float * ffim, * ffim_out;
00048
00049 float
00050 mag,
00051 sum,
00052 sumnz,
00053 mean,
00054 max,
00055 amax,
00056 smax;
00057 int
00058 npts, nznpts;
00059
00060 float wtsum ;
00061 float * wt=NULL ;
00062 PARSER_code * pcode ;
00063
00064 int nw , nnw , iw ;
00065 float * sw=NULL , vw ;
00066
00067 nxy = nx*ny; nxyz = nxy*nz;
00068
00069 #define GOODVOX(ijk) (fmask==NULL || fmask[ijk]!=0)
00070 #define BADVOX(ijk) (fmask!=NULL && fmask[ijk]==0)
00071
00072 ENTRY("EDIT_filter_volume") ;
00073
00074
00075
00076 if( fmask != NULL && filter_opt == FCFLAG_AVER ) filter_opt = FCFLAG_MEAN ;
00077
00078
00079
00080 if( filter_opt == FCFLAG_AVER ){
00081 if( fim_type != MRI_float ){
00082 ffim = (float *) malloc (sizeof(float) * nxyz);
00083 if( ffim == NULL ){
00084 fprintf(stderr,"EDIT_filter_volume: no workspace!\n") ;
00085 EXIT(1) ;
00086 }
00087 EDIT_coerce_type (nxyz, fim_type, vfim, MRI_float, ffim);
00088 } else {
00089 ffim = (float *) vfim ;
00090 }
00091 EDIT_aver_fvol(nx,ny,nz,dx,dy,dz,ffim,filter_rmm) ;
00092 if( ffim != vfim ){
00093 EDIT_coerce_autoscale(nxyz, MRI_float, ffim, fim_type, vfim);
00094 free(ffim) ;
00095 }
00096 EXRETURN ;
00097 }
00098
00099
00100
00101
00102
00103 mask = MCW_build_mask (nx, ny, nz, dx, dy, dz, filter_rmm);
00104 if (mask == NULL)
00105 {
00106 fprintf (stderr, "Warning: Filter option has no effect. \n");
00107 EXRETURN;
00108 }
00109 mnum = mask->num_pt;
00110
00111
00112
00113 if( filter_opt == FCFLAG_EXPR ){
00114 double atoz[26] ;
00115
00116 if( fexpr == NULL ){
00117 fprintf(stderr,"*** EDIT_filter_volume: no fexpr for FCFLAG_EXPR!\n");
00118 EXIT(1) ;
00119 }
00120
00121 pcode = PARSER_generate_code( fexpr ) ;
00122 if( pcode == NULL ){
00123 fprintf(stderr,"*** EDIT_filter_volume: illegal fexpr!\n"); EXIT(1);
00124 }
00125
00126 wt = (float *) malloc(sizeof(float)*(mnum+1)) ;
00127
00128 #define II 8
00129 #define JJ 9
00130 #define KK 10
00131 #define RR 17
00132 #define XX 23
00133 #define YY 24
00134 #define ZZ 25
00135
00136 for( ii=0 ; ii < 26 ; ii++ ) atoz[ii] = 0.0 ;
00137
00138 wt[0] = PARSER_evaluate_one( pcode , atoz ) ;
00139
00140 for (jma = 0; jma < mnum; jma++){
00141 atoz[II] = mask->i[jma] ;
00142 atoz[JJ] = mask->j[jma] ;
00143 atoz[KK] = mask->k[jma] ;
00144 atoz[XX] = atoz[II] * dx ;
00145 atoz[YY] = atoz[JJ] * dy ;
00146 atoz[ZZ] = atoz[KK] * dz ;
00147 atoz[RR] = sqrt(atoz[XX]*atoz[XX] + atoz[YY]*atoz[YY] + atoz[ZZ]*atoz[ZZ]) ;
00148 wt[jma+1] = PARSER_evaluate_one( pcode , atoz ) ;
00149 }
00150 free(pcode) ;
00151 }
00152
00153
00154
00155 if( filter_opt > FCFLAG_WINSOR &&
00156 filter_opt < FCFLAG_WINSOR+FCFLAG_ONE_STEP ){
00157
00158 static int first=1 ;
00159
00160 nw = filter_opt - FCFLAG_WINSOR ;
00161 nnw = mnum - nw ;
00162 filter_opt = FCFLAG_WINSOR ;
00163
00164 fprintf(stderr,"++ Winsor filter: N=%d nw=%d\n",mnum+1,nw) ;
00165 if( first || nnw < nw ){
00166 first = 0 ;
00167 if( nnw < nw ){
00168 fprintf(stderr,"** Illegal Winsor parameters - skipping!\n") ;
00169 EXRETURN ;
00170 }
00171 }
00172
00173 sw = (float *) malloc(sizeof(float)*(mnum+1)) ;
00174
00175 fmask = NULL ;
00176 }
00177
00178
00179
00180 ffim = (float *) malloc (sizeof(float) * nxyz);
00181 if (ffim == NULL)
00182 {
00183 fprintf (stderr, "\n Error: cannot allocate filter workspace! \n");
00184 EXIT(1);
00185 }
00186 ffim_out = (float *) malloc (sizeof(float) * nxyz);
00187 if (ffim_out == NULL)
00188 {
00189 fprintf (stderr, "\n Error: cannot allocate filter workspace! \n");
00190 EXIT(1);
00191 }
00192
00193
00194
00195 EDIT_coerce_type (nxyz, fim_type, vfim, MRI_float, ffim);
00196
00197
00198
00199 for (k = 0; k < nz; k++)
00200 {
00201 for (j = 0; j < ny; j++)
00202 {
00203 for (i = 0; i < nx; i++)
00204 {
00205
00206
00207 ijkvox = THREE_TO_IJK (i, j, k, nx, nxy);
00208 npts = nznpts = 0 ;
00209 sum = sumnz = max = amax = smax = wtsum = 0.0 ;
00210 if( GOODVOX(ijkvox) ){
00211 mag = ffim[ijkvox];
00212 switch (filter_opt)
00213 {
00214 case FCFLAG_MEAN:
00215 sum = mag; npts = 1; break;
00216 case FCFLAG_NZMEAN:
00217 if (mag != 0.0)
00218 {sumnz = mag; nznpts = 1;}
00219 break;
00220 case FCFLAG_MAX:
00221 max = mag; npts = 1 ; break;
00222 case FCFLAG_AMAX:
00223 amax = fabs(mag); npts = 1 ; break;
00224 case FCFLAG_SMAX:
00225 smax = mag; npts = 1 ; break;
00226 case FCFLAG_EXPR:
00227 sum = wt[0]*mag ; wtsum = wt[0] ; npts = 1 ; break ;
00228
00229 case FCFLAG_WINSOR:
00230 for( iw=0 ; iw <= mnum ; iw++ ) sw[iw] = mag ;
00231 vw = mag ; break ;
00232 }
00233 }
00234
00235
00236
00237 switch (filter_opt)
00238 {
00239 case FCFLAG_MEAN:
00240 for (jma = 0; jma < mnum; jma++)
00241 {
00242 ii = i + mask->i[jma];
00243 jj = j + mask->j[jma];
00244 kk = k + mask->k[jma];
00245 if (ii < 0 || jj < 0 || kk < 0 ||
00246 ii >= nx || jj >= ny || kk >= nz) continue;
00247 ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00248 if( BADVOX(ijkma) ) continue ;
00249 mag = ffim[ijkma];
00250 sum += mag; npts++;
00251 }
00252 break;
00253 case FCFLAG_NZMEAN:
00254 for (jma = 0; jma < mnum; jma++)
00255 {
00256 ii = i + mask->i[jma];
00257 jj = j + mask->j[jma];
00258 kk = k + mask->k[jma];
00259 if (ii < 0 || jj < 0 || kk < 0 ||
00260 ii >= nx || jj >= ny || kk >= nz) continue;
00261 ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00262 if( BADVOX(ijkma) ) continue ;
00263 mag = ffim[ijkma];
00264 if (mag != 0.0) {sumnz += mag; nznpts++;}
00265 }
00266 break;
00267 case FCFLAG_EXPR:
00268 for (jma = 0; jma < mnum; jma++)
00269 {
00270 ii = i + mask->i[jma];
00271 jj = j + mask->j[jma];
00272 kk = k + mask->k[jma];
00273 if (ii < 0 || jj < 0 || kk < 0 ||
00274 ii >= nx || jj >= ny || kk >= nz) continue;
00275 ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00276 if( BADVOX(ijkma) ) continue ;
00277 mag = ffim[ijkma];
00278 sum += wt[jma+1]*mag; npts++; wtsum += wt[jma+1] ;
00279 }
00280 break ;
00281 case FCFLAG_MAX:
00282 for (jma = 0; jma < mnum; jma++)
00283 {
00284 ii = i + mask->i[jma];
00285 jj = j + mask->j[jma];
00286 kk = k + mask->k[jma];
00287 if (ii < 0 || jj < 0 || kk < 0 ||
00288 ii >= nx || jj >= ny || kk >= nz) continue;
00289 ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00290 if( BADVOX(ijkma) ) continue ;
00291 mag = ffim[ijkma];
00292 if (npts == 0 || mag > max) max = mag;
00293 npts++ ;
00294 }
00295 break;
00296 case FCFLAG_AMAX:
00297 for (jma = 0; jma < mnum; jma++)
00298 {
00299 ii = i + mask->i[jma];
00300 jj = j + mask->j[jma];
00301 kk = k + mask->k[jma];
00302 if (ii < 0 || jj < 0 || kk < 0 ||
00303 ii >= nx || jj >= ny || kk >= nz) continue;
00304 ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00305 if( BADVOX(ijkma) ) continue ;
00306 mag = ffim[ijkma];
00307 if (npts == 0 || fabs(mag) > amax) amax = fabs(mag);
00308 npts++ ;
00309 }
00310 break;
00311 case FCFLAG_SMAX:
00312 for (jma = 0; jma < mnum; jma++)
00313 {
00314 ii = i + mask->i[jma];
00315 jj = j + mask->j[jma];
00316 kk = k + mask->k[jma];
00317 if (ii < 0 || jj < 0 || kk < 0 ||
00318 ii >= nx || jj >= ny || kk >= nz) continue;
00319 ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00320 if( BADVOX(ijkma) ) continue ;
00321 mag = ffim[ijkma];
00322 if (npts == 0 || fabs(mag) > fabs(smax)) smax = mag;
00323 npts++ ;
00324 }
00325 break;
00326 case FCFLAG_WINSOR:
00327 for (jma = 0; jma < mnum; jma++)
00328 {
00329 ii = i + mask->i[jma];
00330 jj = j + mask->j[jma];
00331 kk = k + mask->k[jma];
00332 if (ii < 0 || jj < 0 || kk < 0 ||
00333 ii >= nx || jj >= ny || kk >= nz) continue;
00334 ijkma = THREE_TO_IJK (ii, jj, kk, nx, nxy);
00335 sw[jma+1] = ffim[ijkma];
00336 }
00337 break;
00338
00339 default: break;
00340 }
00341
00342
00343
00344 switch (filter_opt)
00345 {
00346 case FCFLAG_MEAN:
00347 ffim_out[ijkvox] = (npts > 0) ? sum/npts : 0.0 ; break;
00348
00349 case FCFLAG_NZMEAN:
00350 ffim_out[ijkvox] = (nznpts > 0) ? sumnz/nznpts : 0.0 ; break ;
00351
00352 case FCFLAG_EXPR:
00353 ffim_out[ijkvox] = (wtsum != 0.0) ? sum/wtsum : 0.0 ; break ;
00354
00355 case FCFLAG_MAX: ffim_out[ijkvox] = max ; break;
00356 case FCFLAG_AMAX: ffim_out[ijkvox] = amax; break;
00357 case FCFLAG_SMAX: ffim_out[ijkvox] = smax; break;
00358
00359 case FCFLAG_WINSOR:
00360 qsort_float( mnum+1 , sw ) ;
00361 ffim_out[ijkvox] = (vw < sw[nw]) ? sw[nw]
00362 : (vw > sw[nnw]) ? sw[nnw]
00363 : vw ;
00364 break ;
00365
00366 default: break;
00367 }
00368 }
00369 }
00370 }
00371
00372
00373 EDIT_coerce_autoscale(nxyz, MRI_float, ffim_out, fim_type, vfim);
00374
00375
00376 KILL_CLUSTER (mask);
00377 free (ffim);
00378 free (ffim_out);
00379 if( wt != NULL ) free(wt) ;
00380 if( sw != NULL ) free(sw) ;
00381
00382 EXRETURN;
00383 }
00384
00385
00386
00387
00388
00389
00390
00391 void EDIT_aver_fvol( int nx, int ny, int nz,
00392 float dx, float dy, float dz, float * fim , float rmm )
00393 {
00394 MCW_cluster * mask ;
00395 int i, j, k , ij , ii ;
00396 int jk,jkadd , nxadd,nyadd,nzadd , nxyz_add , mnum ;
00397 float * ffim ;
00398 int * madd ;
00399 float fac , sum ;
00400
00401 ENTRY("EDIT_aver_fvol") ;
00402
00403
00404
00405 mask = MCW_build_mask(nx,ny,nz, dx,dy,dz, rmm) ;
00406 if( mask == NULL || mask->num_pt < 2 ){
00407 fprintf(stderr,"Warning: EDIT_aver_volume has no effect.\n") ;
00408 EXRETURN ;
00409 }
00410 mnum = mask->num_pt ;
00411
00412
00413
00414 #if 1
00415 nxadd = nyadd = nzadd = 1 ;
00416 for( ii=0 ; ii < mnum ; ii++ ){
00417 i = abs((int)mask->i[ii]) ; nxadd = MAX(i,nxadd) ;
00418 j = abs((int)mask->j[ii]) ; nyadd = MAX(j,nyadd) ;
00419 k = abs((int)mask->k[ii]) ; nzadd = MAX(k,nzadd) ;
00420 }
00421 #else
00422 nxadd = (int)(rmm/dx) ;
00423 nyadd = (int)(rmm/dy) ;
00424 nzadd = (int)(rmm/dz) ;
00425 #endif
00426
00427 nxyz_add = (nx+2*nxadd) * (ny+2*nyadd) * (nz+2*nzadd) ;
00428
00429 ffim = (float *) malloc(sizeof(float) * nxyz_add) ;
00430 if(ffim == NULL){
00431 fprintf(stderr,"*** EDIT_aver_volume can't malloc workspace!\n") ;
00432 fprintf(stderr,"nx=%d ny=%d nz=%d nxadd=%d nyadd=%d nzadd=%d\n",
00433 nx,ny,nz , nxadd,nyadd,nzadd ) ;
00434 EXIT(1) ;
00435 }
00436 for( i=0 ; i < nxyz_add ; i++ ) ffim[i] = 0.0 ;
00437
00438 madd = (int *) malloc( sizeof(int) * (mnum+1) ) ;
00439 if( madd == NULL ){
00440 fprintf(stderr,"*** EDIT_aver_volume can't malloc workspace!\n") ;
00441 EXIT(1) ;
00442 }
00443 madd[0] = 0 ;
00444 for( ii=0 ; ii < mnum ; ii++ ){
00445 madd[ii+1] = mask->i[ii] +
00446 mask->j[ii] * (nx+2*nxadd) +
00447 mask->k[ii] * ((nx+2*nxadd)*(ny+2*nyadd)) ;
00448 }
00449 mnum++ ; fac = 1.0 / mnum ;
00450
00451 KILL_CLUSTER(mask) ;
00452
00453
00454
00455 for( k=0 ; k < nz ; k++ ){
00456 for( j=0 ; j < ny ; j++ ){
00457 jkadd = j * (nx+2*nxadd) + k * ((nx+2*nxadd)*(ny+2*nyadd)) ;
00458 jk = j * nx + k * (nx * ny) ;
00459 for( i=0 ; i < nx ; i++ ) ffim[i+jkadd] = fim[i+jk] ;
00460 }
00461 }
00462
00463
00464
00465 for( k=0 ; k < nz ; k++ ){
00466 for( j=0 ; j < ny ; j++ ){
00467 jkadd = j * (nx+2*nxadd) + k * ((nx+2*nxadd)*(ny+2*nyadd)) ;
00468 jk = j * nx + k * (nx * ny) ;
00469 for( i=0 ; i < nx ; i++ ){
00470 sum = 0.0 ; ij = i+jkadd ;
00471 for( ii=0 ; ii < mnum ; ii++ ) sum += ffim[ij+madd[ii]] ;
00472 fim[i+jk] = fac * sum ;
00473 }
00474 }
00475 }
00476
00477 free(ffim); free(madd);
00478 EXRETURN;
00479 }